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Chapter  I 


INTRODUCTION 

Glass  fiber  reinforced  epoxy  is  an  effecient  structural  material, 
coupling  high  strength  with  light  weight.  This  favorable  strength- 
weight  ratio  makes  the  material  attractive  for  some  flight  structures 
as  well  as  other  machines  and  structures  where  weight  is  an  important 
consideration.  In  recent  years  the  material  has  undergone  considerable 
development,  and  it  has  experienced  a  moderate  amount  of  use  for  some 
components.  Presently,  its  use  is  hampered  by  the  difficulty  in 
predicting  the  material's  behavior  under  loads  which  approach  the  breaking 
load.  Whereas  many  materials  exhibit  near-linear  behavior  up  to  failure, 
glass-epoxy  laminates  typically  exhibit  a  considerable  amount  of  nonlinear 
deformation  prior  to  gross  fracture.  Zones  of  increasing  material  damage 
occur  in  the  form  of  crazing,  ply  cracking,  and  ply  delamination .  In 
some  applications,  efficient  use  of  the  material  requires  employing  the 
material  at  stress  levels  well  into  the  nonlinear  portion  of  the  laminate's 
stress-strain  curve.  Questions  then  arise  relative  to  a  precise  definition 
of  exactly  what  constitutes  failure,  the  form  of  the  stress  distribution 
around  notches  and  holes,  and  the  laminate's  stress-strain  response  in  the 
post-crazing  region. 

These  questions  provided  the  impetus  for  this  research.  The  overall 
purpose  of  the  work  was  to  obtain  information  that  would  contribute  to 
rational  methods  of  strength  predictions  and  design  of  glass-epoxies. 

Toward  this  end,  a  joint  program  of  material  testing  and  numerical 


investigation  of  glass-epoxy  laminates  was  undertaken.  The  specific 
objectives  were  to  determine  the  appropriate  stress-strain  behavior 
for  various  glass-epoxy  laminates,  develop  a  finite  element  computer 
model  for  determining  the  stress  distribution  around  stress  raisers, 
such  as  notches  and  holes ,  and  determine  ways  to  apply  these  results 
to  the  failure  mechanisms  to  predict  the  strength  of  laminated  glass- 
epoxy  structures. 

A  specific  material,  Scotchply  XP-250,  was  selected  as  a  test 
material.  This  material  was  characterized  with  respect  to  its  ply 
properties  in  tension,  compression  and  shear.  For  verification  of 
prediction  methods,  three  angle  ply  laminates — [±30]^,  [+45]^,  and 

[±60]  — and  one  quasi-isotropic  laminate,  [0/145/90]  ,  were  tested 

s  s 

in  tension.  Finally,  the  [0/±45/90]s  laminate  containing  a  hole 
was  tested  to  failure.  These  tests  were  compared  to  the  finite 
element  results. 

The  finite  element  program  contains  a  doubly-curved,  thick-shell, 
isoparametric  element.  The  program  will  predict  the  stress  distribution 
in  both  thick  and  thin  plates  and  shells  loaded  transversely  or  inplane. 
The  shape  of  the  crazed  or  yielded  region  around  stress  raisers  can 
be  mapped  and  the  change,  with  increasing  load,  in  the  shape  and  size 
of  this  damage  zone  can  be  determined  on  up  to  complete  failure. 


Chapter  II. 


MATERIAL  CHARACTERIZATION  TESTS 


2 . 1  Introduction 

The  material  chosen  for  this  work  is  known  as  Scotchply  XP-250 
manufactured  by  the  3M  Company.  It  is  a  high-strength,  moldable,  epoxy¬ 
glass  prepreg.  For  the  present  tests  it  was  obtained  in  unidirectional 
cured  sheet  form,  having  either  8  or  14  plys.  The  nominal  ply  thickness 
is  0.009  inch  and  the  fiber  volume  ratio  is  about  50  percent. 

To  obtain  lamina  characterization  of  the  material,  five  types  of 

unidirectional  tests  were  run.  Tension  tests  at  0  and  90  degrees  to  the 

T 

fiber  direction  were  used  to  determine  the  stiffness  properties  E^  , 

T  T  T  T 

E^  ,  and  together  with  the  ultimate  strengths  and  X2  parallel 

and  transverse  to  the  fiber  direction.  Shear  tests  were  used  to  determine 

the  shear  stiffness  G  and  ultimate  shear  strength  Compression 

c  c 

tests  were  used  to  find  the  compressive  stiffness  properties  E^  ,  E^  > 

c  c  c 

and  together  with  the  compressive  ultimate  strengths  X^  and  X^  . 

Specimens  for  all  tests  were  cut  oversize  (1/16  to  1/8  inch)  with 

a  band  saw.  Final  dimensions  were  obtained  by  grinding,  with  water  flowing 

over  the  cutting  area.  All  specimens  were  instrumented  with  350-ohm 

strain  gages.  Because  of  the  small  gage  section,  1/16-inch-long  tee 

rosettes  were  used  on  the  compression  specimens.  Gages  1/4  inch  long  were 

used  on  all  other  specimens.  All  tests  were  conducted  at  room  temperature 

of  about  70°F  and  room  humidity  of  about  50  percent. 


All  specimens  were  tested  in  an  Instron  testing  machine  using  a 
fixed  cross-head  speed.  The  Instron  machine  chart -recorded  load  versus 
time.  Strain  data  were  recorded  by  the  use  of  a  four -channel  strain 
gage  signal  conditioner  together  with  two  two-channel  strip  chart  recorders 
Figure  1  shows  the  test  set  up.  The  load-time  and  strain-time  curves 
are  digitized  by  use  of  a  Tektronics  4051  Computer  Graphics  System.  A 
pen  is  moved  along  the  load-time  curve  taking  load  readings  at  certain 
time  intervals.  The  pen  is  then  moved  along  the  strain-curve,  reading 
strain  values  at  the  same  time  intervals  as  for  the  load  curve.  The 
Tektronics  computer  was  programmed  to  construct  stress-strain  curves  from 
these  data  and  to  compute  the  required  stiffness  parameters. 

2 . 2  Tension  Tests 

Five  tension  tests  were  conducted  on  8 -ply  unidirectional  specimens 

loaded  at  0  degrees  to  the  fiber  direction.  Specimen  dimensions  were 

fixed  by  the  ASTM  standards,  reference  1.  The  specimens  were  9  inches 

long  and  0.5  inch  wide  equipped  with  1. 5- inch- long  load  tabs  made  from 

printed  circuit  board  material.  Tabs  were  attached  before  final  machining 

with  Eastman  910  cement.  Specimens  were  equipped  with  longitudinal  and 

transverse  strain  gages.  Cross-head  speed  was  0.05  inch/minute.  A 

typical  stress-strain  curve  is  shown  in  Figure  2.  The  elastic  modulus 
T 

En  was  determined  from  a  first-order  least  squares  curve  fit  of  only  the 
initial  points  on  the  stress-strain  curve. 

The  failed  0-degree  specimens  are  shown  in  Figure  3  and  the  test 
results  are  summarized  in  Table  1.  The  desired  "shaving  brush"  failure  is 
exhibited  by  specimens  2  and  5  and  by  specimen  3,  although  somewhat 
imperfectly.  Note  from  Table  1  the  resulting  high  ultimate  stress  for 
specimens  2  and  5.  Tab  bond  failure  occurred  on  specimen  1,  possibly 


depressing  its  measured  ultimate  strength.  Ultimate  strength  was  not 

obtained  from  specimen  4  because  a  strain  gage  was  broken  prior  to  the 

ultimate  load  and  the  test  was  stopped  to  investigate.  The  specimen 

was  broken  later.  The  average  values  for  the  elastic  modulus  and 

0 

ultimate  strength  of  5.64  x  10  psi  and  134  ksi  respectively  agree 

0 

well  with  the  values  of  5.70  x  10  psi  and  130  ksi  published  by  the 
3M  Company,  reference  2. 

Five  unidirectional  tensile  tests  were  conducted  at  90  degrees 
to  the  fiber  direction.  The  specimen  dimensions  were  the  same  as  for 
the  0-degree  tests  except  for  the  width  which  was  1.000  inch  instead 
of  0.500  inch,  in  accordance  with  ASTM.  A  typical  stress-strain  curve 
is  shown  in  Figure  4.  The  five  failed  specimens  are  shown  in  Figure  5, 
and  the  test  results  are  summarized  in  Table  2.  The  specimens  6,  7,  and 
10  exhibit  the  desired  type  of  failure,  i.e.  away  from  the  end  tabs. 

The  failure  stress,  however,  of  specimen  9,  which  failed  at  the  tabs,  was 
the  highest  of  all  and  the  failure  stress  of  specimen  11,  which  also  failed 
at  the  tabs,  was  among  the  highest.  This  suggests  that  for  90-degree 
specimens,  failure  near  the  tabs  causes  no  serious  error  in  the  measured 
ultimate  stress. 

T 

Poisson's  ratio  v  was  calculated  from  v,^  =  ^22^11*  ^or 

comparison  the  measured  value  of  is  shown.  The  difference  in  the  two — 

0.078  measured  as  compared  to  0.092  calculated —  is  not  surprising, 

considering  that  the  measured  values  are  obtained  from  a  very  small 

transverse  strain.  Even  a  small  error  in  this  strain  would  account  for 

the  difference  in  the  measured  and  calculated  values.  From  Table  2  it 

T  6 

can  be  seen  that  the  material  has  an  elastic  modulus  of  E ^  =  1.74  x  10 

psi  and  fails  at  an  ultimate  stress  of  7.55  ksi. 


2.3  Shear  Tests 


In  determining  lamina  shear  properties,  the  three-rail  fixture 
was  used.  Shear  testing  has  been  the  subject  of  considerable  controversy 
and  a  number  of  fixtures  or  specimens  have  been  used,  including  the 
+ 45-degree  specimen  (see  references  3  and  4),  the  10-degree  off-axis 
specimen  (see  references  5  and  6),  and  the  two- and  three-rail  fixtures. 
Currently  the  two-  and  three -rail  fixtures  are  being  considered  by  ASTM 
Committee  D-30  as  standard  fixtures  for  finding  inplane  shear  properties. 

Figure  6  shows  two  views  of  the  three-rail  specimen.  The  test 
plate  is  clamped  between  stationary  rails  on  the  edges  while  a  third  rail 
clamping  the  plate  at  the  center,  is  pushed  down  by  the  test  machine, 
loading  the  specimen  in  shear  parallel  to  the  fibers.  Two  strain  gages 
are  attached  to  the  specimen  at  45  degrees  to  the  fiber  direction.  From 
the  strain  transformation  equations  the  normal  strain  at  45  degrees 
is  related  to  the  shearing  strain  referred  to  the  material's  axis  by 

\2l2%5  U) 

The  shearing  stress  between  the  rails  is  assumed  to  be  uniform 

throughout  the  specimen  length  from  top  to  bottom.  So  assumed. 


the  shearing  stress  is  given  by 


12  2bh 

where  P  is  the  load,  b  is  the  specimen  shear  length,  which  was  6  inches, 
and  h  is  the  plate  thickness.  It  is  apparent  that,  while  this  expression 


may  be  accurate,  it  is  not  exact  since  the  shearing  stress  must  go  to 


zero  at  the  top  and  bottom  free  edges  of  the  plate.  Because  of  this, 
the  accuracy  of  the  average  stress,  Equation  (2),  has  been  questioned, 
especially  for  laminates  with  angle  plys.  A  Fourier  series  solution  by 
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Whitney  et  al  [7]  indicated  that  for  a  [1  45]s  laminate  the  shear 
stress  near  the  edge  increased  from  zero  to  a  peak  value  50  percent 
greater  than  the  average  at  a  distance  of  only  one -tenth  the  plate 
length  from  the  edge.  A  recent  finite  element  solution  by  Bergner 
[8],  however,  disagrees  with  this.  This  is  shown  in  Figure  7,  taken 
from  [8].  Bergner' s  results  indicate  that  the  average  shear  distribution 
is  indeed  an  accurate  estimate  of  the  actual  shear  stress  distribution. 

The  present  finite  element  computer  results  (see  Chapter  IV)  are  also 
shown  on  Figure  7.  The  present  finite  elements  results  are  slightly 
lower  than  Bergner' s  but  they  essentially  show  that  a  condition  of 
uniform  shearing  stress  exists  along  the  length  of  the  three-rail  fixture. 

For  the  three-rail  fixture  two  equal  test  sections  exist  on  either 
side  of  the  middle  rail.  Thus  to  fully  utilize  the  specimen,  strain 
gages  are  placed  on  both  sides  and  strain  data  are  recorded  from  both. 

The  three-rail  specimens  must  be  used  with  care.  The  rails  hold  the 
specimen  by  clamping  friction  rather  than  by  bearing  on  the  bolts.  In 
fact,  the  rails  have  emery  cloth  bonded  to  them  and  the  rail  bolts  are 
torqued  to  70  ft-lbs  to  prevent  slipping.  The  holes  in  the  test  plate 
are  considerably  larger  than  the  bolts —  1/2  inch  as  compared  to  3/8 
inch  for  the  bolts.  As  a  consequence,  it  is  possible  to  assemble  the 
specimen  and  fixture  with  considerable  misalignment,  with  the  middle 
rail  tilted  from  vertical,  say.  This  would  destroy  the  assumed  equality 
of  the  test  sections  on  each  side  of  the  middle  rail,  compounded  by  the 
fact  that  the  load  head  would  now  push  down  at  the  top  on  one  side  of 
the  middle  rail  rather  than  at  the  rail's  center.  To  alleviate  this 
problem,  cylindrical  spacers  with  diameters  equal  to  the  width  of  the 
test  section  were  used  to  align  the  rails  during  bolt-up.  These  spacers 


are  visible  in  the  top  picture  of  Figure  6 .  To  further  decrease 
misalignment  the  top  of  the  center  rail  was  machined  so  as  to 
leave  only  a  small  area  1/2  inch  in  diameter  for  the  load  head  to 
push  against. 

Four  specimens  were  tested,  yielding  eight  sets  of  data.  The 

failed  specimens  are  shown  in  Figure  8  and  the  test  results  are 

summarized  in  Table  3.  The  average  lamina  shear  modulus,  was 

0 

0.68  x  10  psi  and  the  ultimate  shear  stress,  ,  was  7.23  ksi.  The 
shear  stress-strain  curves  exhibited  considerable  nonlinearity.  For 
predicting  laminate  behavior  it  was  decided  that  use  would  be  made 
of  the  full  stress-strain  curve  rather  than  just  the  initial  slope.  To 
permit  this ,  all  the  shear  data  were  plotted  and  fitted  with  a  second 
order  least-squares  curve.  Figure  9.  This  curve  now  becomes  the  master 
shear  curve  for  use  in  the  laminate  programs. 


Compression  Tests 

Compression  testing  was  carried  out  using  a  fixture  similar  to  the 
IITRI  compression  fixture,  reference  [9].  Two  views  of  the  fixture 
are  shown  in  Figure  10.  An  exploded  view  drawing  is  shown  in  Figure  11. 
While  the  IITRI  fixture  was  discussed  in  reference  [9]  no  dimensions 
were  given,  and  so  the  fixture  has  been  re-designed  here.  It  was  built 
in-house  requiring  283  man-hours  of  shop  time.  The  fixture  was  built 
from  cold -rolled  steel.  The  two  main  parts  of  the  fixture  are  guided 
together  by  two  rods  0.750  inch  in  diameter  which  fit  into  linear  bearings 
in  the  upper  half.  The  specimen  is  gripped  by  wedges  which  are  bolted 
to  the  specimen  prior  to  the  test.  The  wedge  angle  is  11  degrees.  The 
sloping  surfaces  of  the  wedges  were  lubricated  prior  to  testing  to 
increase  the  wedging  action.  The  wedges  were  2.500  inches  in  length  and 


1.500  inches  in  width.  They  were  slotted  on  the  straight  side  to 
receive  the  tabbed  test  specimen.  These  slots  were  given  gripping 
"teeth"  by  punching  the  slot  surface  numerous  times  with  an  impact 
punch. 

Specimens  were  prepared  for  testing  at  both  0  and  90  degrees  to 
the  fiber  direction.  The  dimensions  of  both  types  of  specimens  were 
the  same.  To  minimize  buckling  the  specimens  were  14  plys  thick.  The 
specimens  were  0.25  inch  wide  and  5.5  inches  long,  equipped  with  crossply 
end  tabs  2.50  inches  long.  This  leaves  a  gage  section  which  measures 
only  0.5  inch  by  0.25  inch — a  small  area  for  a  strain  rosette.  Specimen 
dimensions  are  shown  in  Figure  12.  A  photograph  of  a  specimen  instrumented 
with  strain  rosettes  and  lead  wires  on  both  sides  is  shown  in  Figure  13. 

The  compression  fixture  was  checked  by  conducting  a  test  on  a  2024-T4 
aluminum  specimen.  The  same  specimen  was  then  used  in  a  tension  test 
and  the  stress-strain  curves  for  tension  and  compression  were  compared. 

For  the  compression  test  the  longitudinal  strain  was  monitored  on  both 
sides  of  the  specimen  to  assess  the  degree  of  bending.  The  load-time 
curve  and  the  two  longitudinal  strain-time  curves  are  shown  together  in 
Figure  14.  At  point  A  it  appears  that  minute  grip  adjustment  occurred 
so  that  the  two  strain  curves  abruptly  crossed,  one  increasing  while  the 
other  decreased  with  no  change  in  load.  This  suggests  that  one  of  the 
tabs  on  one  side  of  the  specimen  slipped  slightly  while  its  opposing 
neighbor  held  firm.  This  would  introduce  bending  into  the  specimen 
even  without  fixture  misalignment.  The  strain  on  either  side  deviated 
from  the  average  strain  by  about  6  percent.  The  "flat  spot"  on  the  load¬ 
time  curve  at  B  does  not  indicate  grip  slippage  or  tab  failure  but 
instead  results  from  crosshead  backlash  when  the  testing  machine  is  loaded 
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in  compression. 

The  average  compressive  strain  was  used  to  construct  the  stress- 

strain  curve  shown  in  Figure  15.  That  figure  also  contains  the 

tensile  stress-strain  curve.  Aluminum  is  known  to  possess  very 

nearly  the  same  stress-strain  behavior  in  tension  as  compression • 

Thus  friction  in  the  compression  fixture  would  result  in  a  stress- 

strain  curve  whose  slope  is  too  steep  and  the  measured  value  of  E 

would  be  larger  than  that  for  tension.  The  good  comparison  shown 

in  Figure  15  means  this  does  not  happen,  indicating  that  friction  in 

the  compression  fixture  is  negligible.  The  gripping  problem  reflected 

by  point  A  in  Figure  14,  however,  is  a  source  of  error  which  can 

affect  the  measured  values  of  E  if  the  strain  is  monitored  on  only  one 

side  of  the  specimen.  Moreover,  any  induced  bending  will  tend  to 

depress  the  measured  values  of  the  compressive  ultimate  strength. 

Six  compression  tests  were  run  on  0-degree  specimens.  The 

failed  specimens  are  shown  in  Figure  16.  A  sample  stress-strain 

curve  is  shown  in  Figure  17.  In  that  test — test  17 — as  well  as  two 

other  tests,  strain  gages  were  mounted  on  both  sides  of  the  specimen. 

The  longitudinal  strain  from  both  sides  is  shown  in  Figure  17.  They 

are  not  quite  the  same,  as  ideally  would  be  the  case  if  no  bending 

c 

or  buckling  were  present.  The  two  values  of  E^  obtained  from  both 

6  6 

strains  are  6.39  x  10  psi  and  5.73  x  10  psi.  This  was  the  highest 

difference  in  E^  obtained  from  the  three  specimens  instrumented  on 

both  sides.  A  summary  of  the  0-degree  test  results  is  shown  in 

c  6  . 

Table  4.  For  test  18  the  values  obtained  for  E^  are  5.35  x  10  psi 

6  c  6 

and  5.88  x  10  psi  and  for  test  19  both  values  for  E^  are  5.97  x  10  psi. 

The  strains  from  opposite  sides  of  a  specimen  usually  agreed  fairly  well 


until  some  type  of  failure  (perhaps  fibers  breaking)  began  to  occur. 
This  initial  failure,  marked  by  audible  noise  and  a  small  drop  in  the 
load,  usually  occurred  near  two-thirds  of  the  ultimate  load.  After 
this  initial  damage  or  failure  had  occurred  symmetry  was  lost.  The 
strain  suddenly  increased  in  a  step  fashion  on  one  side  of  the  specimen 
while  suddenly  decreasing  on  the  other  side,  indicating  a  sudden 
application  of  bending  strain.  It  may  be  that  failure  of  a  bundle  of 
fibers  on  one  side  of  the  specimen  causes  a  load  eccentricity  on 
the  remaining  effective  net  section,  hence  bending  must  occur. 

Minute  uneven  slippage  of  the  tabs  in  the  grips  as  already  discussed 
could  cause  the  same  behavior.  Too,  it  must  be  remembered  that  the 
gage  section  is  short  and  that  failure  sometimes  initiated  underneath 
a  gage,  which  could  cause  erratic  gage  behavior.  In  any  case,  the 
data  obtained  after  the  initial  failure — a  sudden  decrease  in  load 
accompanied  by  a  sudden  increase  or  decrease  or  both  in  strain — must 
be  viewed  with  suspicion,  since  the  assumed  symmetry  of  the  test  is 
lost  at  that  point.  For  this  reason  the  measured  ultimate  strains 

Q 

were  not  recorded  in  Table  4 .  The  average  measured  value  of 

G  T 

was  5.87  x  10  psi,  slightly  higher  than  the  value  of  ,  which 

was  5.64  x  106  psi.  Poisson's  ratio,  too,  was  slightly  higher  in 

compression  than  tension — 0.317  as  compared  to  the  tensile  value 

of  0.299. 

Five  compression  tests  were  run  on  90-degree  specimens.  The 
failed  specimens  are  shown  in  Figure  18.  A  sample  stress-strain  curve 
is  shown  in  Figure  19,  and  results  are  summarized  in  Table  5.  More 
than  for  any  other  tests,  the  90-degree  compressive  stress-strain 
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curves  exhibited  an  early  nonlinearity.  Considering  this,  a  secant 
definition  of  might  be  appropriate,  however,  for  the  sake  of 

consistency  with  the  other  tests,  the  slope  of  the  initial  portion 

c  c 

of  the  curve  was  used  for  E^  •  This  made  the  determination  of  E^^ 

0 

difficult,  since  for  this  method,  E ^  depends  strongly  upon  the 

first  few  points  of  the  curve.  This  accounts  for  some  of  the  variation 

in  E^^  shown  in  Table  5.  Bending,  however,  was  a  problem;  test  28, 

monitored  with  gages  on  both  sides,  exhibited  considerable  bending 

0 

as  can  be  noted  from  two  considerably  different  values  of  E  for 
that  test. 
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2 . 5  Summary  of  Material  Properties 

For  handy  reference,  the  various  lamina  properties  determined 
from  the  characterization  tests  for  the  XP-250  material  are 
summarized  below: 


E  T 
L11 

= 

5.64  x  106  psi 

T 

61 

= 

24,000  ue 

« 

E  T 
b22 

= 

1.74  x  106  psi 

x  c 
X1 

= 

112  ksi 

T 

V12 

= 

0.299 

T 

X2 

= 

7.55  ksi 

G12 

= 

0 

0.680  x  10  psi 

T 

S2 

= 

4,760  ue 

• 

E11C 

= 

0 

5.87  x  10  psi 

X  c 
X2 

= 

25.0  ksi 

■ 

E22C 

= 

2.12  x  106 

c 

e2 

= 

18,600  yc 

- 

c 

v ,  „ 

= 

0.317 

S,n 

7.23  ksi 

- 

12 

12 

T 

X1 

= 

134  ksi 

S12 

= 

19,700  ue 

Chapter  III. 


POST-CRAZING  CHARACTERIZATION  OF 
GLASS-EPOXY  LAMINATES 

3.1  Introduction 

A  laminate  contains  a  number  of  laminae  (plys)  oriented  at 
various  angles  to  the  primary  load  direction.  For  loads  limited  to 
the  linear  range,  given  the  lamina  properties  the  usual  lamination 
theory  [10]  gives  accurate  estimates  of  the  overall  stiffness  and 
compliance  of  a  given  laminate.  With  increasing  loads,  however, 
certain  plys  within  the  laminate  begin  to  fail  by  matrix  cracking  and 
splitting  between  the  fibers.  In  glass  epoxies  the  onset  nf  matrix 
cracking  gives  the  laminate  a  hazy,  milky,  light-colored  appearance, 
sometimes  referred  to  as  crazing.  Beyond  the  onset  of  crazing  the 
laminate  compliance  increases  with  increasing  load;  the  crazing  area 
in  a  ply  grows  and  may  extend  to  plys  of  other  angles  before  the 
ultimate  laminate  load  is  reached.  Depending  upon  the  laminate's 
layup,  the  onset  of  crazing  may  occur  at  loads  which  are  rather  low 
compared  to  the  laminate's  ultimate  load.  For  many  structural 
applications  the  laminate's  reserve  strength  beyond  crazing  may  safely 
be  utlized.  For  confident  design  in  much  cases  it  is  important  to  have 
knowledge  of  the  post-crazing  stress-strain  response  of  the  laminate. 

Efforts  at  forming  a  lamination  theory  of  failure  have  genera] ’y 
been  only  moderately  successful.  The  usual  approach  is  to  numerically 
apply  the  laminate  stress  or  strain  in  increments.  After  each  load 
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increment  some  selected  failure  theory  is  applied  to  each  ply.  The 
load  increments  are  continued  until  a  ply  fails,  after  which  the 
stiffness  of  that  ply  is  modified  to  reflect  its  reduced  load 
carrying  capacity.  This  reduced  ply  stiffness  is  then  used  to 
assemble  the  overall  laminate  stiffness  and  the  load  increments 
are  continued  until  other  plys  fail,  after  which  their  stiffnesses 
are  also  reduced.  This  process  is  continued  until,  by  some  definition, 
enough  plys  have  failed  to  constitute  laminate  failure.  This  approach 
has  been  lucidly  discussed  by  Rowlands  [11]  in  the  proceedings  from 
a  ASME  Symposium  on  inelastic  behavior  of  composite  materials  (see 
also  the  Rowlands  report  [12]);  papers  by  Sandhu  [13],  Hahn  and  Tsai  [14] 
and  Chow  et  al  [15,  16]  illustrate  aspects  of  this  approach.  While 
this  approach  is  conceptionally  clear  and  logically  sound,  in  its 
application  a  number  of  problems  must  be  resolved.  In  the  first 
place,  a  ply  may  fail  in  a  number  of  modes — e.g.  splitting  or  crushing 
of  the  matrix  between  the  fibers  due  to  large  transverse  tension  or 
compression,  fiber  failure  in  tension  or  compression,  etc.  How  should 
the  failed  ply's  various  stiffness  constants  be  modified  for  each  mode 
of  failure?  In  other  words,  how  does  the  ply  unload  after  its  failure. 
Evidence  indicates  that  in  situ  ply  strength  and  post-failure  stiffness 
properties  may  vary  considerably  from  those  of  a  unidirectional  test 
coupon  [17],  Futhermore,  a  uniform  definition  of  laminate  failure, 
applicable  to  a  number  of  layups,  is  lacking.  In  some  cases  laminate 
failure  is  assumed  to  occur  once  fiber  failure  (as  distinct  from  matrix 
failure)  has  occurred  in  two  or  more  plys.  This  definition  may  be 
adequate  for,  say,  a  [0/t  45],  layup  but  totally  inappropriate  for  an 
angle  ply  layup  of,  say,  [t  45]s.  other  instants  laminate  failure 


is  assumed  to  occur  once  the  modified  laminate  stiffness  become 


singular.  Each  definition  may  apply,  albeit,  each  to  a  different  class 
of  layups. 

The  above  method  of  laminate  behavior  prediction — here  loosely 
referred  to  as  Rowland's  method  although  a  number  of  researchers  have 
used  it — is  investigated  in  detail  in  the  following.  The  method  is 
assessed  by  applying  it  to  a  number  of  materials — graphite-epoxy, 
boron  epoxy  and  glass-epoxy--having  various  layups.  Biaxial 
failure  response  of  several  glass  epoxies  is  illustrated. 

3.2  Failure  Theories 

A  number  of  failure  theories  are  available  for  predicting  ply 
failure.  An  exhaustive  review  of  failure  theories  for  anisotropic 
materials  was  provided  by  Sandhu  [18];  Rowlands  [11]  also  discussed 
several.  Only  two  were  considered  in  the  present  work:  the  Hill 
theory  [19]  and  the  Tsai-Wu  theory  [20]. 

For  an  orthotropic  ply  in  plane  stress  the  Hill  theory  takes 
the  form. 


°1  °1°2 


(x1)‘ 


12 


=  1 


12 


(3) 


where  and  X^  are  the  uniaxial  strengths  parallel  and  transverse  to 

the  fibers  and  is  the  ply  shear  strengths.  The  Hill  theory  does 

T  T 

not  distinguish  between  the  tensile  strengths  X^  ,  X^  and  the  compressive 

c  c  •  .  «  .  • 

strengths  X^  ,  .  Some  writers  have  made  this  distinction  by  using 

c  T 

X^  when  is  negative  and  X1  when  is  positive  and  similarly  for 
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The  Tsai-Wu  failure  criterion  accounts  for  both  tensile  and 
compressive  strengths.  In  addition  to  quadratic  terms  it  contains 
linear  terms  which  distinguishes  between  negative  and  positive  stresses. 
For  plane  stress  conditions  the  criterion  is  expressed  as. 


Fl°l  +  F2°2  +  F6T12  +  Fll°l  +  F22°2  +  2F 


12°1°2  +  F66T12 


(4) 


where 
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X1  X1 


—  ■  F  = 
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66  S  +S  ' 
12  b12 


1  c 

As  before  ,  X^  are  the  ply  longitudinal  tensile  and  compressive 
T  c 

strengths  and  X2  ,  are  the  ply  transverse  tensile  and  compressive 
strengths.  S12+  and  S10  are  the  ply  strengths  in  positive  and  negative 


12 


inplane  shear.  For  most  composite  materials  S12  =  =  s12  so  that 

Fg  =  0.  The  interaction  term  F^  cannot  be  expressed  in  terms  of  the 

unaxial  strength  properties,  but  must  instead  be  determined  from  biaxial 

tests.  Since  the  accuracy  of  F^2  is  sensitive  to  the  type  of  test  used 

to  find  it,  accurate  values  of  F^2  are  difficult  to  obtain.  Certain 

2  - 

stability  conditions  limit  the  range  of  F^  such  that  F^.j.F22  -  F^2  >  0. 

For  a  glass-epoxy  material  having  the  properties 


X  T  =  154  ksi  X1C  =  88.5  ksi 

*r 

X2  =4.56  ksi  X2C  =  17.04  ksi 


12 


9.00  ksi 


(5) 
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the  limiting  values  of  F  2  become, 

F12  =  1  >/FiiF22  =  ~  9‘717  X  10"10  in4/lb2  (6) 

For  a  unidirectional  off-axis  tension  coupon  the  value  of  F  ^  has  little 
effect  on  the  predicted  strength  as  shown  in  Figure  20.  For  that  case 
the  Tsai-Wu  theory  predicts  about  the  same  strength  whether  F  =  0, 

Ff2  =  +  |/FiiF22  °r  F12  =  ”  ^11^22’  For  reason  the  off-axis  test 

is  known  to  be  an  unsatisfactory  test  for  finding  F^  [21,  22].  Since 
the  influence  of  F^  for  the  off-axis  test  is  small  the  suggestion  is 
that  F  ^  may  be  set  equal  to  zero  without  losing  accuracy.  This  has  been 
common  practice  for  graphite-epoxy  and  will  be  adopted  here  for  glass- 
epoxy  as  well.  Figure  20  also  contains  a  strength  prediction  based 
on  the  Hill  theory.  The  Hill  theory  prediction  agrees  well  the  Tsai-Wu 
prediction  for  the  off-axis  coupon. 

For  angle  ply  test  coupons  the  Hill  and  Tsai-Wu  theories  are 
compared  again  in  Figure  21.  The  value  of  F^  which  gives  the  best 
comparison  iwth  the  Hill  theory  is  F^  =  -  /f'^F^  •  in  the  region, 

0°  >  «  >  35°,  the  Hill  theory  predicts  a  significantly  lower  failure  load 
than  does  the  Tsai-Wu  theory.  Beyond  =  35°  the  two  agree  fairly  well. 
The  value  of  F  2  which  agreed  best  with  the  Hill  prediction  was 
F  =  -  9.717  x  10~10. 

The  Tsai-Wu  failure  theory,  because  of  its  generality  and  because 
it  provides  for  a  difference  in  tensile  and  compressive  strengths  was 
selected  for  the  following  work. 


3.3  Failure  Surfaces  for  Glass-Epoxy  Laminates 


To  illustrate  the  failure  response  to  biaxial  stress,  the  failure 

surfaces  of  several  glass-epoxy  laminates  are  shown  in  Figures  22  thru 

25.  The  failure  strengths  were  predicted  for  various  values  of  the 

stress  ratio  o^/a^  using  the  Tsai-Wu  failure  theory.  The  failure  behavior 

of  a  [±45]^  angle-ply  laminate  is  shown  in  Figure  22.  The  failure 

surface  is  seen  to  be  an  ellipse.  This  is  as  expected  since  the  Tsai-Wu 

failure  theory  applied  to  an  orthotropic  lamina  in  two-dimensional  stress 

space  is  an  ellipse  and  since  the  [+  45]g  laminate  is  essentially  an 

orthotropic  plate.  The  laminate’s  strength  for  hydrostatic  compression, 

the  third  quadrant,  is  great  compared  to  its  strength  for  hydrostatic 

tension,  the  first  quadrant.  In  general  the  laminate's  predicted  strength 

is  great  for  negative  applied  stresses.  The  failure  surface  of  a  [±35]s 

laminate  is  shown  in  Figure  23.  This  laminate  is  stronger  along  the 

x-direction  than  the  y-direction  and  thus  the  long  axis  of  the  failure 

surface  is  skewed  toward  the  axis  in  stress  space.  As  for  the  [±45], 

laminates,  abundant  strength  is  exhibited  in  the  third  quadrant  compared 

to  the  first  quadrant.  Although  complex  structural  shapes  and  complex 

loads  sometimes  result  in  a  laminate  loaded  in  quadrants  2,  3  or  4, 

most  laminates  are  utilized  in  the  first  quadrant  of  the  stress  space. 

That  is,  thin  laminates  are  primarily  tension  structures.  The  failure 

surfaces  for  a  number  of  angle  plys  were  computed  for  the  first 

quadrant  only.  These  are  shown  in  Figure  24.  The  surfaces  for  <*  =  60° 

and  55°  are  the  same  as  for  30°  and  35°  if  a  and  o  are  interchanged. 

x  y 

The  [+15]  and  [+75]  laminates  extend  off  the  graph  exhibiting  considerabl 
s  s 

longitudinal  strength  and  relative  transverse  weakness  on  the  x-  and  y- 


directions  respectively. 


The  character  of  the  failure  behavior  for  a  [0/90 ]  glass -epoxy 
laminate  is  shown  in  Figure  25.  Actually,  three  surfaces  are  shown. 

The  inside  curve — dotted  line — represents  the  first  ply  failure  (FPF) 
in  the  matrix  material  of  one  set  of  plys.  At  this  value  of  the  load 
the  stiffness  properties  of  the  damaged  ply  were  reduced  (see  the  following 
section  for  details)  and  load  application  was  continued  until  matrix  failure 
occurred  in  the  remaining  plys — solid  lines.  The  stiffness  values  of  these 
plys  were  likewise  reduced  and  load  application  was  continued  until 
longitudinal  (fiber)  failure  occurred  in  one  of  the  sets  of  plys;  this 
is  shown  by  the  outside  line.  In  quadrant  one,  great  reserve  strength 
exists  beyond  matrix  failure  of  the  first  two  sets  of  plys.  In  quadrants 
2,  3,  and  4,  however,  the  curve  for  longitudinal  failure  mostly  coincides 
with  the  curve  for  second  matrix  failure;  longitudinal  failure  is 
simultaneous  with  second  matrix  failure. 

3 . 4  Laminate  Response  by  the  Method  of  Rowlands 

Since  the  method  of  Rowlands  [12]  deals  with  stresses  rather  than 
strain  energy  and  uses  the  usual  lamination  theory  the  method  has  strong 
appeal  for  engineers.  The  method  is  concept ionally  simple  and  well  founded 
within  the  framework  of  lamination  theory.  It  was  decided  to  investigate 
this  approach  for  a  wide  range  of  materials  and  layup  configurations. 

The  purpose  was  to  assess  its  applicability  in  general  for  predicting 
laminate  strength  and  to  investigate  its  ability  for  predicting 
specifically  the  strengths  of  glass-epoxy  laminates. 

A  computer  program  was  written  similar  to  that  described  by  Rowland 


in  References  [11]  and  [12].  The  program  predicts  the  inplane  stress-strain 
response  of  a  symmetric  laminate  test  coupon  subjected  to  a  biaxial  test. 
Figure  28.  The  program  contains  both  the  Hill  and  Tsai-Wu  failure  criteria 
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although  only  the  Tsai-Wu  theory  is  used  in  the  following  examples.  The 

laminate  stress  in  the  x-direction  is  applied  in  increments ,  A o  .  The 

laminate  stress  in  the  y-direction  is  given  by  A =  BAo^.  The  average 

shearing  stress  is  given  by  either  At  =  or  by  t  =  constant. 

The  operator  selects  the  values  of  8  and  y  for  a  desired  stress  ratio. 

The  provision  t  =  constant  allows  one  to  obtain  a  failure  curve  for  a 
r  xy 

constant  value  of  shearing  stress . 

For  each  increment  of  stress ,  the  incremental  laminate  strain 

components  Ae^,  and  Ay^  are  calculated  using  the  laminate  compliance 

matrix  from  the  previous  stress  increment.  These  strain  increments  are 

then  used  to  calculate  increments  of  stress  (Ac  ),  ,  (Ao  ),  and  (At  ), 

x  k  y  k  xy  k 

for  each  k  ply  using  the  stiffness  matrices  from  the  previous  load 
increment.  These  stress  increments  are  transformed  to  the  1-2  direction 
for  each  ply  yielding  (Aa2)k,  an<*  ^At12^k‘  current  ply 

stresses  are  then  given  by  adding  the  incremental  stresses  for  the  (n+1) 
cycle  to  the  stresses  for  the  n  load  cycle : 


o^n+1^  =  o^n^  +  Aa^n^ 
o2(n+l)k  =  °2(n)k  +  Aa2(n)k 


Tl2(n+1)k  +  T12^n^  +  Axl2^n\ 


The  total  strains  for  the  n+1  load  increment  are  given  by 


e  (n+1)  =  e  (n)  +  A e 

XXX 


ey(n+l)  =  ey(n)  +  Aey 


y  (n+1)  =  y  (n)  +  Ay 
'xy  xy  xy 


(8) 
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The  average  laminate  stresses,  of  course,  are  given  by 


a  (n+1)  =  a  (n)  +  Acs 

X  X  X 

a  (n+1)  =  a  (n)  +  Acs 

y  y  y 

t  (n+1)  =  t  (n)  +  At 

xy  xy 


xy 


(9) 


The  components  (cs.^,  ^a2^k’  an<*  ^T12^k  are  usec*  in  eit^er  the  Hill 
or  Tsai-Wu  formula  to  investigate  the  failure  of  each  ply.  Once  failure, 
as  predicted  by  the  formula,  is  reached  the  ply  is  next  investigated  to 
determine  if  the  failure  is  matrix  or  fiber  in  nature.  This  is  done  by 
comparing  (cs^)^  with  the  ultimate  tensile  stress  and  ultimate  compressive 
stress.  If  (a  )  exceeds  neither  of  these,  it  is  assumed  that  the  failure 

X  K 

is  in  the  ply  matrix.  After  matrix  failure  if  (a^^  is  positive,  the 
failure  is  designated  as  "RESIN  FAILURE  IN  TENSION."  If  (cs,^)^  negative, 
the  failure  is  designated  as  "RESIN  FAILURE  IN  COMPRESSION."  Once  resin 
failure  occurs,  the  constants  E 22  and  are  set  equal  to  near  zero 

(i.e.,  100)  and  E^  retains  its  original  value.  Actually,  E 22  and  G12  can 
be  modified  differently  for  resin  failure  in  tension  and  resin  failure 
in  compression  although  they  are  modified  the  same  way  in  the  following 
examples.  The  resulting  value  of  v,,^  is  approximately  zero  since 
v0  =  v  E„_/E11.  If  (a..),  exceeds  the  compressive  or  tensile  ultimate 
strength  then  the  failure  is  in  the  fiber  and  it  is  assumed  that  all 
stiffness  of  the  ply  is  lost.  Thus  E.^,  ^22’  anc*  ^12  are  a^  set 
approximately  equal  to  zero.  After  ply  stiffness  is  modified  as  above, 
linear  lamination  theory  is  used  to  calculate  new  values  for  the  laminate 
stiffness  and  compliance  matrices  for  use  in  the  next  load  increment. 
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Once  fiber  failure  has  occurred  in  more  than  one  ply,  then  laminate 
failure  is  assumed  and  computations  are  stopped.  In  addition  to  the 
laminate  stresses  and  strains  the  program  indicates  the  stresses  in 
each  ply,  indicates  the  laminate  load  at  which  a  ply  fails  and  tells 
how  the  ply  failed — i.e.  whether  the  failure  was  a  transverse  failure  in 
the  matrix  or  a  longitudinal  failure  of  the  fibers.  Thus  the  laminate 
stress-strain  curve  is  constructed.  This  curve  is  piece-wise  linear, 
changing  its  slope  at  each  ply  failure.  In  some  laminates  several  plys 
fail  tranversely  so  that  laminate  stiffness  becomes  very  low  and  the 
laminate  compliance  becomes  exceedingly  large,  resulting  in  large 
laminate  strains — strains  of  order  1  or  greater.  This  is  also  taken 
to  be  laminate  failure  since  essentially  all  stiffness  is  lost. 

In  the  following,  the  method  of  Rowlands,  as  explained  above,  is 
compared  with  several  test  results  taken  from  the  literature  for 
graphite-epoxy,  boron-epoxy  and  glass-epoxy  laminates. 

Graphite-Epoxy .  Rowlands  [12]  compared  his  predicted  strain  response 
with  test  results  on  [02/±45]s  graphite-epoxy  loaded  at  several  off- 
axis  angles  to  obtain  various  biaxial  stress  ratios.  As  an  exposition 
of  his  method  as  used  in  the  present  report,  the  experimental  data  from 
five  of  his  figures  are  repeated  here.  Figure  27  shows  the  strain 
response  for  the  0-degree  loading  of  the  [02/t45]s  laminate.  Data  from 
five  tests  are  compared  with  predicted  test  response.  Two  possible  failure 
loads,  76  ksi  and  139  ksi  are  indicated,  one  corresponding  to  transverse 
failure  of  the  +45°  plys  and  one  corresponding  to  longitudinal  failure  of 
the  0°  plys.  Neither  is  very  close  to  the  actual  failure  load  of  100.2  ksi 
(The  present  failure  loads,  76  and  139  ksi,  differ  somewhat  from  those 
stated  by  Rowlands  in  his  report  of  81  ksi  and  184  ksi.  The  reason  why 


is  not  known;  a  difference  in  the  applied  stress  increment  will  cause 
some  variation  but  not  enough  to  explain  the  difference.)  The  point 
to  notice  in  Figure  27  is  that  Rowland's  method  gives  a  good  strain 
response,  but  a  poor  estimate  of  the  laminate's  strength — whether  one 
uses  as  laminate  failure  the  initial  transverse  failure  of  the  145°  plys 
or  the  longitudinal  failure  of  the  0°  plys.  Figure  27  also  shows  the 
Ey  strain  response.  Again,  good  strain  response  is  indicated. 

Response  of  the  [02/±45]s  laminate  loaded  in  the  90-degree  direction 
along  with  the  predicted  response  is  shown  in  Figure  28.  The  responses 
for  a  loading  of  24  degrees  and  45  degrees  respectively  are  shown  in 
Figures  29  and  30. 

Figure  31  shows  a  comparison  of  the  Rowlands  method  with  experimental 

data  obtained  by  Daniel  [6]  for  a  [0/±45/90]  test  coupon.  Only  a  few 
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of  the  computed  points  are  shown.  The  computed  failure  load  of  80  ksi 
compares  fairly  well  with  the  test  value  of  74  ksi.  The  method  slightly 
overestimates  the  stiffness  after  the  computed  transverse  failures  of 
the  90°  and  +45°  plys.  A  definite  change  in  the  slope  of  the  test  curve 
is  easily  seen  in  the  region  where  the  computed  failures  of  the  90°  and 
±45°  plys  occur. 

Figure  32  shows  a  comparison  of  the  Rowlands  method  with  Daniel  data 

for  the  [0/+45/90]  laminate  tested  in  uniaxial  tension  at  30°  to  the 
s 

laminate  axis.  As  computed,  the  75°  and  -60°  plys  fail  nearly  at  the 
same  load — 45  ksi  and  47  ksi,  respectively.  The  method  considerably 
overestimates  the  laminate's  strength;  the  predicted  strength  is  86  ksi 
and  the  tested  strength  was  only  63  ksi. 

Boron-Epoxy.  Data  for  two  boron  epoxy  laminates — [0/90/±45]s  and 
[+45]s  each  tested  at  off-axis  angles  of  15  and  30  degrees — were  taken 
from  Coles  and  Pipes  [23]  for  comparison  with  the  Rowlands  method.  Figure 
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33  shows  the  test  and  computed  results  for  the  [0/90/±45]s  laminate 
loaded  at  15  degrees  off  the  axis.  The  method  grossly  over  estimates 
the  strength.  Test  failure  occurred  at  26  ksi,  before  even  the  first 
computed  transverse  failure  of  any  ply.  Roughly,  the  same  comparison 
is  obtained  for  the  30°  off-axis  test  shown  in  Figure  34. 

The  [*45]  laminate  loaded  15  degrees  off-axis  is  shown  in  Figure 
35.  Test  failure  occurs  there  near  the  load  level  of  transverse  failure 
for  each  of  the  plys,  but  at  a  much  higher  strain — about  20,000  vie  for 
e^  computed  to  4,600  ye  for  a  computed  value  of  e^.  The  computed 
stiffness  is  much  greater  than  the  actual  stiffness.  The  30-degree 
off-axis  case  is  shown  in  Figure  36.  Test  failure  occurs  just  below 
transverse  failure  of  the  plys,  but  again  at  a  higher  strain  than 
indicated  by  the  results. 

Glass-Epoxy.  Test  data  of  Hahn  and  Tsai  [14]  for  a  [0/90^]^  glass- 
epoxy  laminate  and  test  data  for  a  [0/±45/90]s  laminate  from  Chow  et  al 
[15]  were  compared  with  the  method  of  Rowlands.  Figure  37  shows  the 
[0/902]s  laminate.  Good  agreement  is  noted  between  the  test  results  and 
the  predicted  results.  Test  failure  occurs  at  50  ksi  and  the  predicted 
failure  occurs  at  56  ksi.  The  [0/902]  laminate  is  a  particularity 
useful  one  for  studying  the  unloading  behavior  of  transversely  failed 
plys.  In  the  analysis  used  by  Petit  and  Waddoups  [24],  negative  values 
are  assigned  to  certain  stiffness  moduli  of  the  failed  plys,  thus  as  the 
applied  load  is  increased  in  increments  the  failed  plys  gradually  give 
up  their  stress,  redistributing  their  load  to  the  remaining  unfailed  plys. 
The  negative  moduli  values  are  maintained  unti^  the  ply's  stress  approaches 
zero  at  which  time  the  moduli  are  equated  to  zero.  In  the  method  of  Rowland 
as  explained  in  [11]  and  as  used  here,  once  a  ply  fails  by  either  the 
Tsai-Wu  or  Hill  theory  the  transverse  stiffness  moduli  are  equated  to 


zero.  This  means  that  a  failed  ply  does  not  unload  at  all.  As  the 
applied  load  is  increased  in  increments  the  failed  ply  maintains  its 
load  without  either  an  increase  or  decrease  in  the  affected  stress. 

At  the  present,  not  much  is  known  about  the  unloading  response  of  failed 
plys;  the  response  may  depend  upon  the  material  properties,  and  the 
stacking  sequence.  A  number  of  unloading  hypotheses  were  tested  by 
computing  laminate  strain  response  for  various  valued  of  the  moduli 
of  the  failed  plys.  When  E  ,  ^12  an<^  v12  t*ie  'transverse-*-y 

failed  plys  were  equated  respectively  to  E  -0.2  Gi2’  an<^ 

0,  excellent  agreement  resulted.  This  is  shown  on  Figure  38.  While 
the  unloading  factor  of  -0.2  is  only  an  empirical  quantity  having  no 
rational  basis,  the  excellent  agreement  obtained  in  Figure  38  must  be 
regarded  as  a  clue  in  the  unloading  behavior  of  failed  plys.  A  slight 
change  in  slope  of  the  test  curve  can  be  seen  near  the  computed  transverse 
failure  stress  of  the  0-degree  plys. 

Predicted  and  test  results  for  the  [0/145/90]  laminate  are 
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shown  in  Figure  39.  The  prediction  of  either  the  ultimate  strength  or 
the  stiffness  for  the  top  portion  of  the  curve  is  not  as  good  as  for 
the  [0/902]  laminate.  This  is  probably  due  to  the  addition  of  the  145- 
degree  plys  with  their  associated  shearing  stresses.  Lamina  shearing 
stress-strain  response  is  usually  nonlinear. 

Summary  of  the  Comparisons.  The  Rowlands  method  allows  for  nonlinear 
laminate  behavior  by  a  loss  in  stiffness  associated  with  the  transverse 
failure  of  the  various  plys  in  a  laminate.  This  results  in  a  stress- 
strain  curve  piecewise  linear — i.e.  with  a  number  of  slope  changes,  each 
one  corresponding  to  the  failure  of  a  ply.  In  the  comparisons  of  Figures 
27  -39  the  method  is  seen  to  give  reasonable  estimates  of  stiffness 
following  the  initial  ply  failures  for  the  graphite-epoxy  and  glass- 
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epoxy  but  not  boron-epoxy.  Much  less  accurate  is  the  strength  prediction. 

If  one  takes  as  the  definition  of  failure,  longitudinal  fiber  failure  or 
singularity  of  the  laminate  stiffness  matrix  then  the  method  in  general 
considerably  over  estimates  the  failure  strength  of  most  of  the 
laminates.  One  exception  was  the  glass  epoxy  with  the  simple  [O/OO^],, 
layup.  For  angle  plys  loaded  along  the  matrix  axis  the  method  does  not 
permit  sufficient  nonlinearity  because  various  plys  do  not  fail  sequentially 
but  instead  fail  all  at  once  (when  the  +•=  plys  fail  so  do  the  -<*  plys). 

It  was  thought  that  the  laminate  nonlinear  behavior  could  better  be 
accomodated  by  using  the  full  lamina  stress-strain  curve  rather  than 
just  the  initial  slope  of  the  curve.  All  of  the  various  lamina  stress- 
strain  curves  are  reasonably  linear  except  for  the  shear  curve.  Therefore, 
it  was  decided  to  use  the  full  lamina  stress-strain  curve  in  the 
lamination  program.  This  is  discussed  in  the  next  section. 

3 . 5  Laminate  Response  with  a  Nonlinear  Lamina  Shear  Curve 

In  order  to  devise  a  method  of  laminate  strength  and  stiffness 
prediction  which  would  permit  a  higher  degree  of  nonlinearity  it  was 
decided  to  use  the  full  nonlinear  lamina  shear  stress-strain  curve 
together  with  lamination  theory.  This  refinement  was  made  for  the  lamina 
shear  curve  only  since  the  and  strain  responses  are  very  nearly 
linear  and  the  shear  strain  response  is  usually  highly  nonlinear. 

The  lamina  shear  curve  was  used  in  a  manner  similar  to  that  of 
Sandhu  [13].  The  actual  shear  stress-strain  curve  was  approximated  by 
a  cubic  spline  function.  This  function  was  incorporated  into  the  program 
described  in  Article  3.4.  Using  the  full  curve,  after  each  increment 
of  stress  the  laminate's  compliance  and  stiffness  are  evaluated  by 
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using  the  tangent  modulus  corresponding  to  the  current  slope  of  the 
ply  shear  stress-strain  curve.  Thus,  general  nonlinear  laminate  response 
is  allowed  over  the  full  range  of  applied  load  values,  including  the 
region  prior  to  first  ply  failure.  The  resulting  laminate  response, 
still  exhibiting  sharp  changes  in  slope  at  each  ply  failure,  will  now 
be  nonlinear  between  the  neccessive  ply  failures  and  not  piecewise 
linear  as  before.  Except  for  making  use  of  the  full  shear  stress-strain 
curve  the  method  is  the  same  as  explained  in  Article  3.4. 

3.6  Test  Laminate  Response  Compared  with  Predicted  Response  Using 
Ply  Nonlinear  Shear  Behavior 

Three  angle  ply  laminates  and  one  quasi-isotropic  laminate  of 

Scotch  Ply  XP-250  were  tested  to  failure  in  uniaxial  tension.  The 

laminate  layups  were,  [130]  ,  [^45]  ,  [t60]  ,  and  [0/±45/90]  .  Three 

s  s  s  s 

tests  were  run  for  each  layup.  The  coupon  dimensions  were  the  same  as 
those  used  for  the  90°  unidirectional  material  characterization  tests, 

1  inch  wide  by  9  inches  long  with  end  tabs  for  gripping.  Strain  gages 
were  used  to  record  the  longitudinal  and  transverse  strains,  and  . 

The  stress-strain  response  of  each  layup  was  determined  and  compared 
with  the  response  predicted  by  the  method  explained  in  Article  3.5. 

The  tensile  stiffness  properties  of  Article  2.5  were  used  in  the 
predictions. 

Figure  41  shows  the  stress-strain  response  of  the  [130]  ^aminate. 

The  response  is  nonlinear  almost  from  the  beginning.  The  predicted 
longitudinal  strain  is  somewhat  greater  than  the  measured  strain  although 
the  difference  is  generally  less  than  10  percent.  The  agreement  for  the 
transverse  strain  is  not  as  good.  As  a  result  of  using  the  ply  nonlinear 
stress-strain  curve  the  computed  curves  in  Figure  41  exhibit  correctly 


the  decreasing  stiffness  with  increasing  load.  For  example,  the  initial 

0 

stiffness  E  of  the  laminate  is  about  3.34  x  10  psi,  whereas  the 
xx  ^ 

0 

stiffness  at  the  predicted  failure  load  is  only  about  2.02  x  10  psi. 

The  predicted  ultimate  stress  is  low,  about  42  ksi  as  compared  to  an 

actual  failure  stress  of  about  60  ksi.  Generally,  strength  predicitions 

using  lamination  theory  fall  below  the  actual  strength  for  angle  plys . 

Chamis  and  Sullivan  [17]  have  indicated  that  this  may  be  due  to  the 

difference  in  the  in  situ  ply  strength  and  the  ply  strength  measured  in 

unidirectional  coupons.  Use  of  the  ply  nonlinear  shear  curve  improves 

the  failure  prediction  only  slightly.  The  failed  [±30]s  coupons  are 

shown  in  Figure  42.  As  can  be  seen,  final  fracture  resulted  from  a 

combination  of  matrix  splitting  between  fibers,  fiber  fracture,  and 

delamination.  Delamination,  indicated  by  the  light  region  around  the 

fracture  surface,  was  extensive.  Final  failure  was  sudden,  with  complete 

loss  of  load  occurring  almost  instantaneously. 

Figure  43  shows  the  initial  portion  of  the  stress-strain  curve 

for  the  [i45]s  laminate;  the  full  curve  is  shown  in  Figure  44.  As 

noted  by  Rotem  and  Haskin  [25]  the  [+45]s  laminates  exhibits  a 

singular  amount  of  large  deformation  prior  to  ultimate  failure.  In 

Figure  44  it  can  be  seen  that  this  laminate  yields  at  a  stress  of  about 

17  ksi.  The  specimen  deforms  by  a  scissoring  action  and  the  strain 

continues  with  little  increase  in  load  to  a  strain  of  about  35,000- 

40,000  ye.  Then  the  curve  starts  climbing  again  and  failure  finally 

occurs  at  a  strain  of  near  100,000  ye — a  10  percent  elongation--four 

times  the  failure  strain  of  the  [160]  laminate.  The  transverse  strain 
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was  practically  equal  to  the  longitudinal  strain  after  the  onset  of 


extensive  yielding.  This  transverse  contraction  is  visible  in  the 

pictures  of  Figure  45  by  comparing  the  width  of  the  tested  specimen 

with  the  width  of  the  end  tabs,  which  originally  were  the  same  width 

as  the  specimen.  As  the  scissoring  action  took  place  crazing  spread 

over  the  whole  specimen  (still  in  evidence  by  the  light  appearance  of 

the  specimens).  Failure  occurred  by  a  combination  of  delamination, 

fiber  breakage  and  splitting  between  the  fibers.  Failure  occurred 

near  the  end  tabs  where  the  scissoring  action  was  restrained  by  the 

stiffness  of  the  tabs.  Figure  43  shows  the  predicted  strain  response. 

The  transverse  and  longitudinal  strains  both  agree  well  with  the  test 

values  up  to  the  predicted  failure  load  of  12  ksi.  The  tangent  modulus 

0 

of  the  curve  decreases  from  about  2.07  x  10  psi  at  the  origin  to 
about  1.13  x  106  psi  at  the  predicted  failure  load  of  12  ksi.  The 
predicted  failure  load  is  too  low,  and  the  extensive  straining  beyond 
17  ksi  followed  by  a  rising  curve  is  not  predicted.  While  the  extensive 
strain  ability  of  the  [+45]s  laminate  is  interesting,  for  most  structural 
applications  the  laminate  could  not  be  utilized  beyond  the  17  ksi  knee 
because  of  the  large  deformations  and  material  damage  associated  with  a 
higher  stress.  It  is  felt  that  from  a  structural  viewpoint  the  useful 
strength  of  the  laminate  is  about  17  ksi  rather  than  the  higher  figure. 

The  predicted  and  test  response  of  the  [l60]s  laminate  is  shown 
in  Figure  46.  While  the  correct  trend  is  predicted,  the  overall  predicted 
stiffness  of  the  laminate  is  greater  than  the  test  stiffness.  Transverse 
failure  of  all  plys  occurs  at  a  stress  of  about  9  ksi.  The  actual 
failure  stress  was  about  11  ksi.  The  predicted  tangent  modulus  decreased 
from  an  initial  value  of  1.70  x  10^  to  a  final  value  of  1.47  x  10^  psi. 
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The  fractured  specimens  are  shown  in  Figure  47.  Practically  no 
delamination  occurred  on  these  specimens  with  the  exception  of  a 
narrow  region  adjacent  to  the  fracture  surface.  Rather,  the  fracture 
extends  along  the  fibers  of  one  set  of  the  plys ,  breaking  the  fibers 
of  the  other  set.  Failure  occurred  by  matrix  splitting  in,  say, 
the  +60-degree  plys  and  by  fiber  failure  in  the  -60-degree  plys. 

The  strain  response  of  the  [0/±45/90]s  laminates  is  shown  in 

Figure  48.  Transverse  failure  of  the  90-degree  plys  is  predicted  at 

a  stress  of  14  ksi  followed  by  a  transverse  failure  of  the  +45-degree 

plys  at  a  stress  of  18  ksi.  Final  failure  is  predicted  when  longitudinal 

fiber  failure  occurs  in  the  0-degree  plys  at  a  laminate  stress  of  53  ksi. 

The  actual  test  failure  stress  was  about  41  ksi.  In  contrast  to  the 

case  of  the  angle  plys,  for  the  [0/145/90]  laminate,  the  prediction 

method  over  estimates  the  strength.  From  a  design  viewpoint  the  method 

erred  on  the  side  of  safety  for  angle  plys  but  for  the  [0/t45/90]s 

laminate  the  results  are  nonconservative.  The  stiffness  of  the  laminate 

is  predicted  very  well,  however.  The  predicted  longitudinal  stiffness 

0 

decreases  from  an  initial  predicted  value  of  3.01  x  10  psi  to  a  final 
value  of  1.88  x  .  06  psi.  The  strain  response  was  also  computed  by 
setting  E22>  G^*  and  v2i  a^ter  t^'e  ply  failure  equal  to  -0.2  times 
their  original  values.  The  prediction  is  shown  in  Figure  49.  The 
agreement  on  the  failure  load  is  improved  to  a  value  of  50  ksi.  The 
stiffness  of  the  upper  portion  of  the  curve  seems  slightly  low,  however. 

It  has  been  noted  before  that  the  ply  unloading  factor  of  -0.2  resulted  in 
good  comparisons  for  some  other  materials.  The  failed  coupons  are  shown 
in  Figure  50.  Half  of  the  45-degree  plys  failed  by  matrix  splitting,  the 
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other  half  by  fiber  failure.  The  90-degrees  plys  of  course  failed  by 
matrix  splitting  and  the  O-degree  plys  failed  by  longitudinal  fiber 
failure.  A  considerable  amount  of  delamination  can  be  seen.  Crazing 
due  to  transverse  failure  of  the  matrix  of  the  90-degree  plys  can  be 
seen  throughout  the  length  of  the  coupon. 

3 . 7  Conclusions  on  Glass-Epoxy  Laminate  Response 

In  the  present  lamination  method,  after  matrix  failure  the  constants 
^22’  ^12  w*iere  set  approximately  equal  to  zero  under  incremental  loading. 
Using  negative  values  for  these  constants  (simulating  unloading)  resulted 
in  no  substantial  improvement  of  predicted  and  test  laminate  stress-strain. 
Use  of  the  nonlinear  ply  shear  curve  resulted  in  a  better  laminate  stress- 
strain  curve — one  which  was  nonlinear  between  failures  of  the  various  plys 
and  also  nonlinear  prior  to  first-ply-failure.  Ply  failure  stresses  and 
strains  agreed  well  with  abrupt  changes  in  the  slopes  of  the  glass-epoxy 
test  curves. 

Two  definitions  of  laminate  failure  were  used:  (1)  longitudinal  fiber 
failure  in  two  or  more  plys  (2)  strains  of  order  one  (singularity  of 
stiffness  matrix).  Definition  (1)  is  appropriate  only  if  a  high  percentage 
of  the  fibers  correspond  to  the  load  direction.  This  definition  over 
estimates  the  ultimate  load  in  the  [0/t45/90]s  laminate  by  about  30  percent. 
For  angle  plys  loaded  along  the  principal  axis,  definition  (1)  does  not 
apply.  All  plys  fail  in  the  matrix  simultaneously  leading  to  very  high 
strains  (singular  laminate  stiffness  matrix)  of  100  percent  or  more  on 
the  very  next  load  increment.  Hence,  definition  (2)  was  used.  This 
resulted  in  failure  predictions  for  angle  plys  which  were  low  by  20  to  30 
percent . 


Improvements  in  the  prediction  of  laminate  ultimate  loads  are 
desirable.  Lamination  theory  is  inherently  limited,  omitting  inter¬ 
laminar  shear  behavior  and  making  no  distinction  in  stacking  sequence. 
Given  these,  it  may  be  that  no  rational  refinement  will  result  in  any 
further  improvement  in  lamination  prediction  of  ultimate  loads.  The 
use  of  in  situ  ply  strengths  as  suggested  by  Chamis  and  Sullivan  [17] 
may  in  the  future  be  a  fruitful  approach. 

For  the  time  being,  the  present  method — using  the  nonlinear  ply 
shear  curve  together  with  E 22  and  G^2  equal  near  zero  after  matrix 
failure — yields  a  laminate  stress-strain  curve  sufficiently  accurate 
for  many  practical  engineering  applications.  This  stress-strain  behavior 
will  be  employed  in  the  finite  element  program  and  example  problems  of 
Chapter  V. 


Chapter  IV. 


DEVELOPMENT  OF  THE  FAILURE  ANALYSIS  METHOD— 

A  DOUBLY-CURVED,  ISOPARAMETRIC,  THICK-SHELL  FINITE  ELEMENT 

4.1  Introduction 

Early  theory  on  laminated  plates  and  shells  [26]  was  a  direct  exten¬ 
sion  of  the  classical  thin  plate  and  shell  theory  based  on  the  so-called 
Kirchhoff  assumptions.  Later  the  bending-extension  coupling  was  studied 
by  Reissner  and  Stavsky  [27].  In  1971  Pryor  and  Barker  [28]  developed  a 
rectangular  finite  element  for  laminated  plates.  The  shear  deformation 
was  included  by  the  relaxation  of  part  of  Kirchhoff' s  assumptions.  Later, 
in  1976,  a  quadrilateral  element  for  laminated  plates  was  presented  by 
Nopratvarakorn  [29].  The  latter  element  is  similar  to  but  more  versatile 
than  the  one  developed  by  Pryor  and  Barker.  The  plate  quadrilateral  ele¬ 
ment  is  then  further  extended  to  model  the  shell  structure.  A  plate  ele¬ 
ment  to  model  shells  has  the  merit  of  simplicity,  but  a  large  number  of 
elements  are  needed  for  modeling  shell  structures.  Therefore,  a  doubly- 
curved,  isoparametric,  quadratic,  8-node,  thick-shell  element  is  developed 
in  this  study.  The  element  is  derived  from  the  16-node  solid  element  by 
specializing  the  element  so  that  strain  energy  of  the  stresses  normal  to 
the  midsurface  is  ignored  and  by  constraining  lines  initially  normal  to 
the  midsurface  to  remain  straight.  Thus  fewer  degrees  of  freedom  are 
needed  to  define  the  displacement  field.  The  resulting  element  has  40 
degrees  of  freedom — three  displacements  and  two  rotations  for  each  of  the 
eight  nodes.  Though  the  midsurface  normals  are  to  remain  straight  during 
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deformation,  these  lines  need  not  remain  normal  to  the  deformed  mid¬ 
surface.  Therefore,  the  ability  to  model  transverse  shear  deformation 
is  retained.  Transverse  shear  is  thought  to  be  significant  for  the 
laminated  plates  and  shells. 


4.2  Isoparametric  Elements 

Considering  the  geometry  of  the  three-dimensional  element  in  Figure 
51,  one  notes  that  by  means  of  the  coordinate  transformation 


x  =  T,  N'x 
i  i  i 

y  =  ?N!y. 

ill 

z  =  EN’z 
i  i  i 


(10) 


the  element  can  have  curved  boundaries.  This  is  an  important  advantage 
of  the  isoparametric  formulation.  In  Equation  (10)  x,  y,  and  z  are  the 
coordinates  at  any  point  of  the  element  and  x^,  y^,  z^,  i  =  1,  .  .  .  n 
are  the  coordinates  of  the  n  nodes.  The  interpolation  functions  N.  are 
defined  in  the  natural  coordinate  system  of  the  element,  which  are  func¬ 
tions  of  5,  n»  C  that  each  vary  from  -1  to  +1. 

In  the  isoparametric  formulation  the  element  displacements  are 
interpolated  in  the  same  way  as  the  geometry;  i.e.,  one  assumes 


u  =  EN'u 
i  i  i 

v  =  EN'v 
i  i  i 

w  =  ZN'w  (ID 

i  i  i 

where  u,  v,  and  w  are  the  local  element  displacements  at  any  point  of 
the  element  and  u..,  v^,  and  w^,  i  =  1,  .  .  .  n,  are  the  corresponding 
element  displacements  at  its  nodes. 

For  a  16-node  solid  element  the  interpolation  functions  are 


r 


r 


r 


r 


r 


f 


t 


I 


defined  to  be 
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N'  =  l/8(l-C)(.l-n)a-c)(-e-n-l) 
N’  =  l/4(l-£2)CL-n)CL-C) 

Ng  =  l/8(l+5)(l-n)(l-?)(5-n-l) 

=  1/4(1- n2)(.l+C)d-c) 

=  i/8(i+c)(i+n)(i-c)(  C+n-D 
Ng  =  i/4(i-c2)(i+n)d-c) 

N}  =  l/8(l-s)(l+n)(l-0(-S+n-i) 
Ng  =  i/4(i-n2)(i-C)(i-c) 


N^  through  N|g  can  be  obtained  by  replacing  5  with  -£.  With  the  defini¬ 
tion  of  N!,  the  first  of  Equations  (10)  can  be  written  as 

8  ( l-r)  ®  (l+t) 

x  =  E  N . ■■■%;-■  x .  +  Z  N.  ■  X4' x.  Q 
i=l  1  2  1  i=i  1  2  1+8 


x  =  i!1Mi(¥L)xiP  +  J1Ni(iri)xi< 


where 


N.  =  i/4(i+eei)(i+nni)(«i+nni-l)  for  i  =  1,  3,  5,  7 
N:  =  l/2(l-C2)(l+nni)  for  i  =  2,  6 
N.  =  l/2(l+«.)(l-n2)  for  i  =  4,  8 

are  the  shape  functionsof  the  8-node  two-dimensional  element  of  the  mid¬ 
surface  ,  and 


ei  =  -1,  0,  1,  1,  1,  0,  -1,  -1 
=  -1,  -1,  -1,  0,  1,  1,  1,  0 


for  i  =  1,  2,  ...  8 


and  x.  ,  x.  ,  etc.,  are  global  cartesian  coordinates  of  the  16  nodes  on 
lq  ip 

z,  =  -1  and  +1.  Similar  expressions  can  be  written  for  y  and  z;  i.e.. 


x.  I  x. 

+  ih"i(¥>kp 


Following  Ahmad  [30]  the  full  three-dimensional  element  is  then  reduced 
to  the  conventional  representation  by  midsurface  nodes  only,  preserving 
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fc 


most  of  the  desirable  characteristics  of  the  solid  element.  The  six 
degrees  of  freedom  can  be  transformed  into  three  mid-point  translations 
and  two  mid-point  rotations  about  two  axes  perpendicular  to  the  normal , 
and  a  change  of  length  of  the  normal  itself.  This  yields  a  stiffness 


too  high  in  bending  due  to  the  fact  that  the  normal  strain  e  =0.  How- 

n 


ever,  Ahmad  replaced  the  linear  £  variation  of  the  normal  displacement 


with  the  condition  =  0,  the  usual  assumption  for  beam  and  plate  theory. 


A  linear  assumption  in  the  £  direction  for  the  in-plane  displacements  u 
and  v  is  sufficiently  good  to  represent  membrane  strain  states  exactly 
and  transverse  shear  strains  closely;  therefore,  all  desired  features 
are  now  included.  Introducing  the  following 


X. 

1 

X.  +x. 

_  ip  iq 

2 

x .  -x . 

y. 

y.  +y. 
ip  J  xq 

V.. 

iq  ip 

y-  -y.  V 

i 

2 

3i 

■’'iq  Jip 

z. 

z.  +z. 

_  ip  iq 

z.  -z. 

1 

2 

iq  ig 

into  Equation  (13)  yields  a  set  of  equations  which  define  element  geometry 


in  terms  of  midsurface  nodal  coordinates  and  vectors  , 

O  1 


v =  iiiw 


‘1  e 
. >  +  .e.n.c* 

if  1=1  1 7rV„. 


(14) 


z . 
L  i 


2  3i 


The  global  coordinate  system  x,  y,  z  has  z vertically  upwards.  The 
local  5,  n>  C  system  is  defined  by  the  intrinsic  shell  coordinates.  The 
five  degrees  of  freedom  of  each  node  will  be  three  translations  u,  v,  w 
in  the  global  x,  y,  z  system  and  two  rotations  a,  6  about  axes  ^ . 
The  directions  of  $  ^  are  t*ie  tangential  direction  of  the  £,  n  coor¬ 
dinates  respectively.  is  drawn  in  the  C  direction  of  the  c-axis  to 
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form  a  triad  with  and  V2  at  the  node,  Figure  52. 


Consider  a  rotation  a.  about  axis  V„.  and  B.  about  axis  V  ..  The 

i  2i  i  li 

displacement  at  any  point  C  from  the  midsurface  is 
?ti,0iVli  6iV2i, 


where  t.^  is  the  thickness  of  the  shell  at  node  i.  Therefore  the  complete 
40-degree-of -freedom  element  displacement  field  may  be  written  as 


hh 


1 =  iIAVi  * 


where 


l  -£ 
li  2i 


=  »u  -m2i  =  »i21  Ui22 


"li  _n2i 


Uill  Pil2 


Ui31  Pi33 


and  ^2i’m2i’n2i  are  resPectively  the  direction  cosines  of 

V1  and  V ^  at  node  i. 

The  strain-displacement  transformation  may  be  obtained  by  differen¬ 
tiating  Equation  (15).  Since  the  displacement  field  is  defined  in  the 
local  Sn?  system,  the  derivatives  in  this  system  must  be  evaluated  before¬ 
hand.  Now  at  any  point 

3u  )y  |w  3^  5z  9u  9v  3w 

ac  as  9S  3?  as  3C  3x  ax  ax 

8u  3v.  3w  =  _8x  3z.  3u  3v  3w 

an  an  an  3n  3n  3n  3y  3y  3y 

3u  3v  3w  3x  3y  3z  3u  3v.  8w 

3C  3c  3c  3  c  3s  3c  3"z  3z  3z  (16) 


or  symbolically 


[uvw.  ,]  =  [J][uvw  ] 
SnS  xyz 


where  [J]  is  the  Jacobean  of  the  transformation  of  x,  y,  z  to  S,  n,  C» 
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Thus 


^xyz3  =  f'ETHuvn^: 


(17) 


Now  the  strains  e..  are  various  combinations  of  » 

13  dx  dX  oX  dy  dy 

3w  3u  3v  3w  ,  ,  .  ,  ,  _  ,  r  t 

g^r,  g-p  -g^-,  —  and  can  be  picked  out  of  the  matrix  [uvw^  zJ  term-by-term 


or 


fe  1 

X 

1 

0 

0 

0 

0 

0 

0 

0 

e 

y 

0 

0 

0 

0 

1 

0 

0 

0 

e 

<  Z 

►  = 

0 

0 

0 

0 

0 

0 

0 

0 

Yxy 

0 

1 

0 

1 

0 

0 

0 

0 

Yyz 

0 

0 

0 

0 

0 

1 

0 

1 

.V 

_0 

0 

1 

0 

0 

0 

1 

0 

w. 


w. 


’y 


W, 


(18) 


Substituting  Equations  (16),  (17)  into  Equation  (18)  leads  to 

{e}  =  CB]{«}  (19) 

where  {6}  is  the  nodal  displacement  matrix.  Matrix  [B]  is  a  6  x  40  and 


[V  = 


of  ei 

ght 

6x5 

blocks.  A  typical  block  [B.]  is  of  the  form 

a. 

1 

0 

0 

t. 

Aiuill-2 

t. 

t . 

Aipil2~T 

t . 

0 

b. 

1 

0 

BitJi21_2 

t. 

Vi22~? 

t. 

- 

• 

0 

0 

c . 

1 

Cipi31~2 
t . 

t. 

,Vi2l'T) 

CiPi32~T 

t  .  t  . 

b. 

1 

a . 

1 

0 

(Vin-^ 

^  Bi*Jil2~2+AilJi2?-'2^ 

’  • 

0 

c . 

1 

b. 

1 

(CiMi21-2H 

(CiMill_2H 

t . 

'ViaiT* 

t . 

(Ciyi22'TfBiUi3?~2) 
t.  t. 

c . 

1 

0 

a . 

1 

ykiviZl^ 

(Ciyil2-2+Aipi32~) 

(20) 

and  [J*]  =  [J] 


A. 

1 


V  +  J$3  "i 


I 


B. 

1 


V 


+ 


N. 

1 


C. 

1 


+  J33 


N. 

l 


4.3  The  Elasticity  Matrix 

Consider  each  lamina  or  each  layer  of  the  composite  behaving  as  a 
homogeneous  orthotropic  material.  Nine  independent  elastic  constants 
are  required  to  describe  the  material.  For  the  principal  axes  of  elastic 
symmetry  (010203),  Figure  53,  which  coincide  with  the  reference  axis 
(x'y'z1)*  the  compliance  relations  for  a  typical  layer  of  a  composite 


are 


r  \ 

(9) 

-  1 

el 

E1 

c* 

^21 

E2 

"  E1 

V31 

e3 

< 

>  = 

'  E1 

Y12 

0 

Y23 

0 

Y!3 

0 

0 

1 

J23 

0 


0 

1 

G!3 


r - 

t— i 

_ 

0 

2 

a 

3 

< 

> 

T 

12 

T 

23 

T 

13 

(0) 


(21 ) 


where  the  coefficient  matrix  is  symmetric.  If  the  fiber  arrangement 
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were  such  that  the  variation  of  the  properties  in  the  2-3  plane  were  to 
be  negligible,  the  number  of  constants  reduces  to  five. 


ei  =  el 


E2  =  E3  '=  ET 


V21  =  V31  VTL 


V  s  V 

23  TT 


G12  =  G13  GLT 


G23  =  GT  =  2(1+vtt) 


(22) 


E  and  G  are  the  Young's  modulus  and  shear  modulus, 
is  defined  as 

e . 

v.  .  = - — 

13  £j 


Poison's  ratio  v. . 

i] 


due  to  a  stress  in  the  j -direction.  Therefore,  the  determination  of  the 
five  independent  elastic  constants  EL,  E^,,  GLT,  and  vTT  character¬ 

izes  the  transversely  isotropic  composite.  If  the  expression  for  G^, 
does  not  hold,  the  number  of  independent  elastic  constants  increases  to 
six. 

The  constitutive  relationship  for  the  generally  orthotropic  com- 

0 

posite  is  obtained  by  solying  for  (o  }  in  Equation  (21)  and  making  use 
of  Equations  (22). 


where  is  neglected  as  in  classical  lamination  theory 


C.,  -  if  i,  i  =  l,  2 

Q.  .  =  3  33 

i] 


C.  . 
i: 


if  i,  j  =  4,  5,  6 


For  an  arbitrary  orientation  of  a  lamina,  the  principal  axes  of 
material  will  not  coincide  with  the  reference  axes  of  the  laminate. 
The  transformation  for  expressing  stresses  in  an  (910263)  coordinate 
system  in  terms  of  stresses  in  (x'y’z’)  system.  Figure  54,  is 


fv 

(e) 

r  9 

COS^(j> 

sin2<|> 

2sin<t>cos<f> 

0 

0 

- 

\o  ,  1 

X* 

°2 

sin2(() 

cos2<|> 

-2sin4>cos<j> 

0 

0 

0  . 
y' 

I12 

>  = 

-sin<]>cos<j>  sin<|>cos$ 

cos2i(i-sin2!t) 

0 

0 

< 

T  ,  , 
x’y 

T23 

0 

0 

0 

cos  <t> 

sin 

<p 

Ty’  z* 

X13 
V-  J 

(K) 

0 

0 

0 

-sin  <t> 

cos 

$ 

Tx'  z’ 

^  J 

(K) 


or 


{0  }(K)  =  [Ta]{a' }(K) 


Hence,  solving  Equation  (26)  yields 


where  [T  J  1  is  obtained  by  replacing  sin<J>  with  -sin$  in  [T  ]. 
o  o 


Similar  transformations  for  the  strain  can  be  obtained  as 
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sin2$  sin4>cos4> 

cos2(J>  -sin$cos<J> 

2sin4>cos<f>  cos2(j!-sin2<J> 
0  0 

0  0 


or 

(£\k>  =  tT£3(E'1(K) 

and 

{e’>(K)  = 

where  [T^]  ^  is  obtained  by  replacing  sin<j>  with  -sinijj  in  [T^]  and  the 
prime  represents  the  principal  axes  system. 

Equations  (25)  can  now  be  rewritten  as 

{°  }(K)  =  [Q](K){e9}(K)  (28) 

Substituting  Equations  (26),  (27)  into  (28)  yields 

<°U(K)  =  CT£:TCQ](K)[Teac)(K)  (29) 

or 

{0*}(K)  =  CQ,;!(K){e,}(K) 

where 

CQ':(K)  =  [T£]TM3(K)[Te3 

and 

T  -1 

[T  J  =  [T  ] 
e  a 


b 
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Equation  (.29)  into  Equation  (30),  the  element  strain  energy  can  be 
obtained  as 


N 


U*.E%J  {e,}T[Q']^,{e’}dV 


K=1 


(K) ' 


(K) 


(31) 


(K) 


where  { e  * }  is  equal  to  since  the  distribution  of  strain  is 

assumed  to  be  continuous  throughout  the  entire  thickness  of  the  lamin¬ 
ate.  Matrix  [Q']^^  is  the  stress-strain  matrix  for  the  layer.  One 
must  ensure  that  [Q]^^  provides  for  zero  stress  normal  to  the  shell. 

Let  the  reference  shell  coordinate  x'y'z'  have  the  same  directions  as 

“k  — ►  — ► 

Vf,  V2,  Vg  so  that  at  each  point  of  the  shell  z'  is  normal  to  the  mid¬ 
surface.  Taking  for  example  the  Kth  layer,  the  stress-strain  relation 


This  form  of  [Q']/„\  provides  for  a  ,  =  0  and  plane  stress  conditions  in 
the  x'y'  plane.  A  coordinate  transformation  is  applied  to  convert 
{e'}(K)  to  matrix  {e}^  in  the  global  xyz  coordinates;  i.e.,  substituting 
{€’}  =  [T»J it} 


into  Equation  (31)  yields 


u  =  Vi  /,  ^)T^pT[QC(K)cr:(E 


}dV 


(K) 


(33) 


(K) 

Introducing  Equation  (19)  into  (33)  leads  to 
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u  =  biarckju} 


(34) 


where 


N 


[k]  =  ^  /  CBJ  [EJ(K)[B3dV(K) 

V(K) 


(35) 


is  the  element  stiffness  matrix  and 


[E]  =  [T'^CQ'D^xCT^ 


(K)1 


in  which  [T*]  is  the  transformation  matrix  between  xyz  and  x’y'z'  coordi¬ 


nates  . 

The  integral  of  Equation  (33)  is  evaluated  by  numerical  integration 
with  respect  to  the  local  n,  5  coordinates.  Matrix  [B]  may  be  split 
into  a  part  [B^]  independent  of  £  and  a  part  D  linear  in  The 
products 


CE](K)CV 


and 


c[B1]  Ce](k)Cbq] 

are  linear  in  C,  representing  the  bending-membrane  coupling  effect.  The 
product 


CV  “VW 


and 


CV  [E]00CB1] 


are  the  membrane  and  bending  effects  respectively. 


4.5  Body  Loads,  Surface  Loads 

Nodal  loads  resulting  from  body  force  and  surface  pressure  will  now 
be  considered.  The  nodal  loads  associated  with  these  applied  forces  may 
be  found  by  usual  procedures  and  only  an  outline  is  included  here . 


Equation  (15)  can  be  rewritten  as  follows: 
{f}  =  jvj  =  tN]{6} 


(36) 


J 


..  'I 


1 
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where  [N]  defines  the  nature  of  the  displacement  field.  Matrix  [N  ]  is 
a  3  x  40  and  is  built  of  eight  3x5  blocks.  A  typical  block  [  N^]  is  of 


the  form 


N.  0  0 

l 


[N.]  =  0  N.  0 

l  l 


0  N. 

l 


nXt  £..  N.£t.£_. 

l  l  li  l  l  2i 

2  "2 

N.£t  .m.  .  N  Xt  .ra0 . 

l  l  li  l  l  2i 

2  ”2 

N.-^t.n,  .  N.£t.n„. 

l  l  li  l  l  2i 

2  2 


The  array  of  element  nodal  forces  {r}  produced  by  body  force  and  surface 

pressure  in  the  element  is  [31J 

{r}  =  /  [N]T{P}dV  +  / [N]T{p  }ds  (37) 

Vol  S  n 

The  first  integral  represents  the  body  force  and  the  last  integral  repre¬ 
sents  surface  normal  pressur  . 


4.6  Computer  Implementation 

First  of  all,  consider  the  definition  of  the  three  mutually  perpen¬ 
dicular  vectors  ,  V^,  ^3i  33  shown  in  Figure  52.  The  rotation  vectors 
0.,  ou  are  colinear  with  v  ^  and  respectively.  It  is  conceivable  that 
in  the  assembled  structure  no  two  nodal  rotation  vectors  will  have  the 
same  direction.  Vector  .  may  be  defined  by  input  data,  and  is  presumed 

O  1 

to  span  the  thickness  and  be  normal  to  the  midsurface.  This  proves  to  be 
very  time-consuming  in  preparation  of  data  for  a  large-scale  problem.  In 
this  study  the  following  approach  is  adopted. 

From  the  differential  geometry,  the  tangent  vector  e^,  e2  as  shown  in 
Figure  55  along  the  local  intrinsic  shell  coordinates  axes  can  be  found 
by  the  following  equations. 


->■  _  9x  -t  ^  3y  -t  3z  > 
e,  =  -7T=-  l  +  1  +  rr  K 


'1  H 


H  J  H 


-*•  3x  f  ^  3y  +  .  3z  i 
e2  "  3n  an  3  3n  K 


•  • 


* 


The  direction  cosines  of  e.  give  the  directions  for  tf...  \T_  .  could  be 

3  3i  li 

defined  by  e^  or  might  be  defined  by  input  data  so  that  it  coincides  with 

a  principal  direction  of  an  orthotropic  material  and 
^ 

V  .  =  V  *  V 
2i  3i  li 

Next  the  element  stiffness  matrix  given  by  Equation  (35)  will  be 
integrated  numerically  with  respect  to  the  local  £,  n,  5  coordinates 


resulting  in 


.1  .1  .1 


M  =  Z  j  f  /  [B]  [E3, „,mdet|j|dCdndc 


-1  -1  -1 


When  the  Jacobian  of  the  above  equation  is  computed,  it  is  found  that  5 
to  the  first  power  appears  in  certain  terms.  These  terms  may  be  neglected 
in  comparison  with  terms  to  which  they  are  added.  In  this  study,  these 
terms  have  been  temporarily  suppressed  so  that  [J]  becomes  independent 
of  t  and  explicit  integration  through  the  thickness  is  possible;  and,  as 
indicated  in  the  previous  section,  [B]  may  be  split  into  [BQ]  +  cCB^]. 

This  integral  results  in 

M  =  f1  /  (CB03Ce3CboMb03Cde3[B1>Cb1][de][bo] 

+[B1J[dJ[B13)x det j J jd£dn 

where  [E]  is  the  in-plane  stiffness  array,  [D]  is  the  flexural  rigidity 
array,  and  [DE]  is  the  coupling  of  membrane  and  bending  stiffness  array. 

The  average  values  of  [E],  [DE],  and  ID]  are  then  computed  and  the  C 
terms  are  restored  in  [J].  Two  Gaussian  points  in  the  thickness  direc¬ 
tion  are  used  for  numerical  integration;  i.e. 

[k]  =  J1  J1  /1([B0][E][Bo>[Bo][DE][B1]+[b1][DE][Bo] 


+[B13[D][B1]  X  det|j|d5dndc 


■  * 
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The  cost  of  numerical  integration  is  then  doubled  but  the  stiffness 
obtained  is  much  better  than  the  one  obtained  by  neglecting  t,  terms  in 


In  the  final  calculation  of  element  stresses,  the  strains  computed 


{e}  =  [B]{6> 

are  referred  to  global  coordinates  x,  y,  and  z.  The  strain  is  then  trans¬ 
ferred  back  to  the  x'jy'jz’  coordinates  by 

{£’}  =  [T']{e} 

And,  finally,  the  operation 

{a'}(K)  =  [Q,:](K){e'} 

gives  stresses  referred  to  the  shell  coordinates  x^y'jZ1;  and  the  stresses 
in  the  principal  material  direction  are  obtained  by 

<Ak)  = 

4.7  Yield  Criteria 

The  Hill  criterion  as  well  as  the  Tsi-Wu  tensor  criterion  are  imple¬ 
mented  in  the  computer  program  to  assess  the  effect  of  stresses  and  strains 
on  the  structural  integrity  of  the  composite.  Due  to  the  thinness  of  the 
plates  and  shells  in  the  93-direction,  plane  stress  is  assumed;  that  is 

°3  ~  T13  ~  t23  s  0 

Also  for  a  composite  which  has  a  regular  fiber  array  in  the  G^G2  plane 


it  is  usually  assumed 


(o„)  =  (o0) 

3  y,p,  2  y.p. 


With  these  conditions  the  Hill  and  Tsi-Wu  criteria  are  essentially  the 
same  as  used  in  the  plane  stress  analysis.  Incorporation  into  the  com¬ 
puter  code  makes  it  possible  to  assess  the  structural  integrity  of  the 
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plates  and  shells  at  discrete  points  due  to  the  finite  element  approxi¬ 
mation.  At  these  discrete  points  the  layer  stresses  a^,  t  £»  T13* 

^22  are  computed  as  described  previously,  and  are  substituted  into  either 
Equation  (3)  or  Equation  (4).  The  layer  is  yielded  if  the  left-hand  side 
of  the  equation  is  greater  than  or  equal  to  1. 

4.8  Mesh  Generation 

The  preparation  of  element  data  is  a  very  time-consuming  task. 
Incorrect  element  data  is  also  a  major  source  of  errors  when  running 
finite  element  programs.  The  mesh  generation  subroutine  is  developed 
to  generate  the  element  data  automatically.  The  MESH8  subroutine  uses 
a  group  of  either  8-node  (quadratic)  or  12- node  (cubic)  quadrilateral 
regions  to  define  the  body  under  consideration.  This  sub-program  is  cap¬ 
able  of  modeling  two-  or  three-dimensional  plates  and  shells  midsurface 
domains  that  are  composed  of  8-node  quadrilateral  elements.  The  element 
nodes  are  numbered  and  the  element  nodal  connectivities  also  generated. 
The  8-node  quadrilateral  region  is  available  in  MESH8.  It  can  be  used  to 
generate  a  two-  or  three-dimensional  quadrilateral  element  with  eight 
nodes.  The  eight  nodes  that  define  the  region  are  numbered  as  shown  in 
Figure  56.  Node  1  is  always  at  the  coordinate  location  £  =  n  =  -1. 

The  region  is  then  subdivided  into  elements  by  considering  eight 
nodes  that  form  a  quadrilateral  such  as  the  area  in  Figure  56  with  the 
center  node  being  omitted. 

The  size  of  the  elements  can  be  varied  by  placing  nodes  2,  4,  6,  or 
8  at  some  point  other  than  the  center  of  the  side .  Movement  of  these 
nodes  shifts  the  origin  of  the  £-n  coordinate  system. 


A  domain  is  generally  modeled  by  using  several  regions  connected  to 
one  another  along  one  or  more  sides.  The  possibility  of  a  common  bound¬ 
ary  between  two  regions  requires  that  proper  connectivity  data  be  given 
These  connectivity  data  convey  to  the  computer  how  the  region  under  con¬ 
sideration  is  connected  to  other  regions. 


Chapter  V 


LAMINATE  STRESS  ANALYSIS  BY  THE  FINITE  ELEMENT  MODEL 

5.1  Description  of  the  Computer  Program 

The  element  is  a  general  doubly-curved  8  node  laminated  thick-shell 
isoparametric  element  which  can  be  used  to  model  both  thick  and  thin 
plates  and  shells,  laminated  or  single  layer.  The  material  can  be 
homogeneous  isotropic  or  orthotropic.  Both  geometric  and  material  non¬ 
linearity  have  been  considered.  The  incremental  procedure  is  employed. 
The  load  increments  are  of  equal  magnitude.  The  load  is  applied  one 
increment  at  a  time  and  during  the  application  of  each  load  increment  the 
equations  are  assumed  to  be  linear.  The  coordinates  of  the  node  are  then 
updated  and  the  adjusted  coordinates  are  used  in  the  computation  of  the 
stiffness  for  the  next  increment.  The  shear  stress-strain  curve  of  the 
composite  material  is  highly  nonlinear.  The  current  tangent  modulus  is 
used  in  the  calculations.  The  shear  curve  is  fit  by  a  cubic  spline 
interpolation  to  define  the  shear  modulus  at  a  given  strain.  Figure  57 
is  a  flow  chart  showing  the  sequence  of  operations  performed  by  the 
program.  Appendix  A  together  with  Figures  58  and  59  give  an  explanation 
of  the  data  input  for  the  program.  Appendix  B  gives  the  program  listing. 

5 . 2  Verification  of  the  Computer  Model 

This  section  presents  the  solutions  of  several  problems  which  are 
intended  to  illustrate  the  capabilities  and  limitations  of  the  finite 
element  computer  program.  Examples  included  have  known  solutions,  and 
thus  provide  good  test  cases  for  the  program. 
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Homogeneous  Simply  Supported  Square  Plate.  The  first  example  is 
the  well  known  Reissner  thick  plate  problem.  A  nondimensional  deflection 
parameter  was  calculated  for  various  plate-thickness-to-lateral-dimension 
ratios.  These  results  are  presented  in  tabular  form  in  Table  6  and  in 
Figure  60.  They  show  excellent  agreement  with  the  Reissner  theory 
(see  [29]  or  [32]).  As  the  thickness  to  lateral  dimension  ratio  gradually 
increases,  the  solutions  of  both  the  finite  element  method  and  Reissner 
theory  disagree  with  the  classical  solution. 

Cylindrical  Shell  Roof.  This  is  a  test  example  of  application  of 
the  element  to  a  shell  in  which  bending  action  is  severe,  due  to  supports 
restraining  deflection  at  the  ends.  The  shell  is  supported  on  diaphragms 
as  shown  in  Figure  61.  These  allow  no  displacements  in  their  own  plane, 
but  offer  no  resistance  to  displacements  perpendicular  to  it.  Only  a 
quarter  of  the  shell  was  actually  analyzed,  by  using  symmetric  boundary 
conditions  along  the  two  orthogonal  planes  of  symmetry.  Displacements  of 
the  shell  in  the  vertical  direction  at  the  mid-span  section  are  shown  in 
Figure  62.  The  reference  curve  is  that  used  by  Pawsley  [33].  The  graphs 
show  that  this  shell  roof  is  well  modelled  by  even  one  element. 

Thin  Hyperbolic  Paraboloid  Shell.  The  boundary  of  this  shell  is 
assumed  to  be  rigidly  held  against  both  displacements  and  rotations. 

The  shell  is  subjected  to  a  uniform  load.  The  geometry  and  material 
properties  are  shown  in  Figure  63.  The  entire  hyperbolic  paraboloid 
was  modelled  using  only  4  elements.  The  results  obtained  are  presented 
along  with  the  results  of  Minch  and  Chamis  [34].  These  results  show 
good  agreement. 

These  comparisons  all  indicate  a  high  degree  of  accuracy  for  the 
present  method.  The  method  will  now  be  applied  to  a  glass-epoxy  laminate 


with  a  hole. 
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5 . 3  Response  of  a  [0/t45/90]s  Glass-Epoxy  Laminate  with  a  Hole 

An  example  like  that  of  Chow  et  al  [16]  was  chosen.  Three  tensile  T" 

coupons  of  XP-250  containing  a  hole  were  tested.  The  layup  was  [0/±45/90]^, 
eight  plys  thick.  During  the  load  application  the  strain  was  monitored 
near  the  hole  by  a  1/16-inch  strain  gage.  The  coupon  dimensions  and  gage  I' 

location  are  shown  in  Figure  65. 

Figure  66  shows  the  mesh  layout  for  the  computer  simulation  of  this 
problem.  The  properties  used  in  the  input  are  those  of  Article  2.5.  The  §'' 

tensile  values  of  the  stiffness  properties  were  used  in  the  input  together 
with  the  nonlinear  ply  shear  curve,  Figure  9.  After  matrix  failure  in  a 
given  ply  E ^  and  G ^  were  set  approximately  equal  to  zero  as  explained  • 

in  Article  3.7.  The  computed  response  is  compared  with  the  three  test 
responses  in  Figure  67.  The  two  agree  fairly  well  although  test  strain 
is  slightly  larger  than  the  computed  strain.  The  indicated  computed  i 

failure  was  taken  to  be  when  two  plys  failed  by  fiber  fracture.  In  this 
problem  these  failures  occurred  first,  of  course,  in  the  elements  on  the 
hole  edge.  The  slight  disagreement  in  Figure  67  is  probably  due  more  ® 

to  imperfect  material  characterization,  as  discussed  in  Article  3.7, 
than  to  the  numerical  method.  The  present  finite  element  model,  as 

already  seen,  appears  to  be  quite  precise.  * 
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Summary  of  Ply  Properties  from  [0°]s  Tension  Tests  on  XP-250  Glass -Epoxy 


Table  3 


Summary  of  Ply  Shear  Properties  for  XP-250 
Glass-Epoxy  from  Three-Rail  Shear 
of  Unidirectional  Panels 


Properties 

Test  No.  's. 

Shear  Modulus 

Gj2  (psi) 

Ultimate  Shear 
Stress 

S12  (ksi) 

Ultimate  Shear 
Strain 

e12  (m  strain) 

12  (I) 

0.60  x  106 

6.A8 

18,200 

12  (II) 

0.63  x  106 

1A, 300 

13  (I) 

0.69  x  106 

7.60 

20,300 

13  (II) 

0.76  x  106 

18,200 

1A  (I) 

0.58  x  106 

6.76 

19,200 

1A  (II) 

0.78  x  106 

15,000 

15  (I) 

0.7A  x  106 

8.  10 

25,900 

15  (II) 

0.66  x  106  i 

26,300 

A 

0.68  x  106 


Mean  Values 


7.23 


19,700 


Table  4 


Summary  of  Ply  Compressive  Properties  for  XP-250 
Glass-Epoxy  from  Tests  of  [0°]  Coupons 


vssProperties 

Test  No. 

Young' s  Modulus 

(psi) 

Poisson's  Ratio 

c 

V12 

Ultimate  Stress 

X*  (ksi) 

16 

4.94  x  106 

0.403 

_ 1 

17 

(Channels  1,  2) 

6.39  x  106 

0.294 

122 

17 

(Channels  3,  4) 

5.73  x  106 

0.309 

18 

(Channels  1,  2) 

5.35  x  106 

0.277 

1 

18 

(Channels  3,  4) 

5.88  x  106 

0.325 

19 

(Channels  1,  2) 

5.97  x  106 

0.330 

121 

19 

(Channels  3,  4) 

5.97  x  106 

0.287 

20 

5.19  x  106 

0.324 

99.8 

21 

6.38  x  106 

0.305 

105 

Mean  Values 

5.87  x  106 

0.317 

112 

^No  ultimate  load  due  to  grip  slippage. 
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Typical  Stress-Strain  Response  for  [0JS  Tensile  Specimen 
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The  Shear  Stress-Strain  Response 


Figure  10.  Two  Views  of  the  Compression 
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Figure  12.  Dimensions  of  Compression  Specimen 


Figure  15.  The  Measured  Compressive  Strain  Response  Compared  with  the  Tensile  Strain 
Response  for  Aluminum  2024-T4 
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igurel7.  A  Sample  Stress-Strain  Curve  for  a  [0]s  Specimen  Under  Compressive  Loading 


Figure  18.  The  Failed  [90]s  Compressive  Spec 


OFF-AXIS  ANGLE, 0-DEGREES 


Figure  20.  Predicted  Strength  of  Glass-Epoxy  Off-Axis  Unidirectional 
Tension  Coupon  by  the  Tsai-Wu  and  Hill  Failure  Theories 
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t igure  21.  Predicted  Strength  of  Glass- Epoxy  Angle  Ply  Tension  Coupon 
by  Tsai-Wu  and  Hill  Failure  Theories 
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Figure  22.  Failure  Envelope  for  a  ( ± 45 ] s  Class-Epoxy  Laminate 


Failure  Envelope  for  a  [+35] s  Glass- Epoxy  Laminate 


Figure  25.  Failure  Envelope  for  a  [0/90] s  Glass- Epoxy  Laminate 
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Figure  27.  Strain  Response  for  the  O-degree  Loading  of  a  [02/±45]s  Graphite- Epoxy  Laminate 


e  for  the  90-degree  Loading  of  a  (02 / 5 ] s  Graphite- Epoxy 
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Figure  29.  Strain  Response  for  the  24-degree  Loading  of  a  [ 0 2 /  ^ 5  J 
Graphite-Epoxy  Laminate 
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Figure  31.  Strain  Response  for  the  0-degree  Loading  of  a  [0/ ±45/901 s 
Graphite -Epoxy  Laminate 
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Figure  33.  Strain  Response  for  a  [0/90/±45]s  Boron-Epoxy  Laminate  Loaded  at  15  Degrees 
Off-Axis 


»ure  36.  Strain  Response  for  a  [ ± A3 ] s  Boron- Epoxy  Laminate  Loaded  30  Degrees  Off-Axis 
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Figure  37.  Strain  Response  of  a  [ 0/ 90 2 1 s  Class- Epoxy  Laminate 
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Figure  38.  Effect  of  Failed-Ply  Unloading  on  the  Computed  Response  of 
the  (0/902 1 s  Glass -Epoxy  Laminate 


STRESS 


Figure  39.  Predicted  and  Test  Strain  Response  for  a  fO/ *  45/901 s 
Glass-Epoxy  Laminate 
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Figure  40.  Effect  of  Failed-Ply  Unloading  on  the  Computed  Response  of 
the  [0/±45/90)s  Glass-Epoxy  Laminate 
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Figure  42.  The  [+30]s  XP-250  Glass- 
Coupons  After  Failure. 
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Figure  47.  The  [l!60]  XP-250  Glass-Epoxy  Tens 

Coupons  After  Failure 
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Figure  49.  Effect  of  Failed-Ply  Unloading  on  the  Computed  Strain  Response  of  the 
[0/±45/90]s  XP-250  Glass~Epoxy  Laminate 


Figure  50.  The  [0/±45/90]s  XP-250  Glass-Epoxy  Tens 
Coupons  After  Failure 


Figure  53.  The  Lamina  Principal  Axes  of  Elastic  Symmetry 


Figure  54.  Relationship  of  Lamina  Principal  Axes 
to  Shell  Coordinates 


Figure  55.  Intrisic  Shell  Coordinate  Axes 


Figure  56.  Possible  Region  for  the  Quadrilateral 
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Figure  58.  The  Sides  of  a  Quadrilateral  Region 


Figure  59.  A  Connected  Set  of  Quadrilateral  Regions 
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Figure  60.  Influence  of  Transverse  Shear  on  Maximum  Deflection  of  a 
Homogeneous  Simply  Supported  Plate 
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Figure  64.  Vertical  Deflection  Across  the  Midspan  of  a  Clamped 
Hyperbolic  Paraboloid  Under  Uniform  Load 
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Figure  65.  Coupon  Dimensions  for  the  [0/+45/90]  Class-Epoxy  Laminate 
with  a  Hole  s 


Figure  66.  Mesh  Layout  for  the  [ 0/^45/90 ]  Coupon  with  a  Hole 


Figure  67.  Comparison  of  the  Computed  and  Test  Strain  Near  a  Hole  in 
a  [0/145/90]  Glass-Epoxy  Laminate 


Appendix  A 


Data  Input  for  the  Program 


Card  1 


Card 


Card 


TITLE  (13A6 ) 

Alphanumeric  statement 

"MESH  ONLY"  if  only  the  mesh  generation  is  desired 


NNP,  Number  of  nodal  points 
NEL,  Number  of  elements 

NMAT,  Number  of  different  filamentary  composite  material 
NSHELL,  0  for  shell  analysis,  1  for  plate  analysis 
NFLAG,  1  for  Hill  criterion,  2  for  Tsai-Wu  criterion 


LD,  Load  identification  1  for  uniform  load  in  x,  y,  z, 
direction,  2  for  uniform  normal  pressure,  3  for  non- 
uniform  normal  load,  4  for  concentrated  load,  5  for 
edge  pressure 

NELPL,  Number  of  element  with  uniform  or  non-uniform  load 
NEDGEL,  Number  of  edge  pressure  boundary  conditions 
NLOAD ,  Number  of  concentrated  loads 
INCMAX,  Maximum  number  of  increments 


Col. 

1-70 

Col. 

70-78 

2  General  1 

Col. 

1-5 

Col. 

6-10 

Col. 

11-15 

Col. 

16-20 

Col. 

21-25 

Leave 

Col.  . 

3  Load  Con' 

Col. 

1-5 

Col. 

6-10 

Col. 

11-15 

Col. 

15-20 

Col. 

20-25 

If  INCMAX  is  set  equal  to  1  linear  analysis  will  be  executed. 


Card  4  Mesh  general  control  card  (215) 

Col.  1-5  INRG,  Number  of  region 

Col.  6-10  INBP,  Number  of  boundary  points  to  be  input 
Card  5  X-coordinates  of  the  boundary  nodes  (8FI0.0) 

Card  6  Y-coordinates  of  the  boundary  nodes  (8FI0.0) 

Card  7  Z-coordinates  of  the  boundary  nodes  (8FI0.0) 


This  format  is  repeated  until  all  the  nodal  values  are  read. 


Card  8  Regions  connectivity  data  (515) 

Col.  1-5  NRG,  Region  number 

Col.  6-10  Four  connectivity  numbers  for  a  region  one  for  each  side. 
Col.  11-15  Each  value  is  the  number  of  the  region  connected  to  a 
Col.  16-20  particular  side.  The  sides  of  the  quadrilatural  region 
Col.  21-25  are  labeled  as  shown  in  Figure  58. 

See  example  for  the  determination  of  the  connectivity  data. 
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Card  9  Region  data  (1115) 

Col.  1-5  NRG,  Region  number 
Col.  6-10  NROWS,  Number  of  rows  of  nodes 
Col.  11-15  NCOL,  Number  of  columns  of  nodes 
Col.  16-20 
Col.  21-25 

Col.  26-30  NDN  Global  node  numbers  used  to  define  the  quadrilateral. 

Col.  31-35 

Col.  36-40 

Col.  41-45 

Col.  46-50 

Col.  51-55 

Replace  Card  4  through  Card  9  with  the  following  data  Cards  if  mesh 
generation  is  not  used.  (Substitute  cards  are  indicated  by  an  asterick.) 

Card  4  Nodal  coordinates  card  (Al,  14,  5X,3F10.0,  15)  (one  for  each  node) 

Col.  2-5  N,  Node  number 

Col.  6-10  Leave  blank 

Col.  11-20  X(N)  X-coordinate 

Col.  21-30  Y(N)  Y-coordinate 

Col.  31-40  Z(N)  Z-coordinate 

Col.  41-45  KN  Node  number  increment 

Nodal  coordinate  card  need  not  be  input  in  node  order  sequence, 
however,  all  nodal  coordinates  must  be  defined.  Joint  data  for  a  series 
of  nodes  may  be  generated  from  information  given  on  two  cards  in  sequence: 


Card 

1 

Nl . 

Card 

CM 

N2 . 

CM 

KN2  is  the  mesh  generation  parameter  given  on  the  second  Card  of 
the  sequence.  The  first  generated  node  is  N^  +  (lx  KN2);  the  second 
generated  node  is  Nq  +  (2  x  KN2);  etc.  Generation  continues  until  node 
number  N2  -  KN2  is  established.  Note  that  the  node  difference  N2  -  Nj_ 
must  be  evenly  divisibly  by  KNj. 

Card  5  Nodal  connectivity  card  (915)  (one  for  each  element) 


Col. 

1-5 

N,  Element 

number 

Col. 

6-10 

NOD 

(N,l) 

Col. 

11-15 

NOD 

(N,2) 

Col. 

16-20 

NOD 

(N,3) 

Col. 

21-25 

NOD 

(N,4) 

Global  nodal  point  number  corresponding 

Col. 

26-30 

NOD 

(N,5) 

to  element  nodes 

Col. 

31-35 

NOD 

(N,6 ) 

Col. 

36-40 

NOD 

(N,7 ) 

Col. 

41-45 

NOD 

(N,8) 

I 


r 


r 


r 
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Card  10  Nodal  ID  Card  (715)  (one  for  each  node) 

Col.  1-5  N,  Node  number 

Col.  6-10  ID(N,1)  x-translation  boundary  condition  code 

Col.  11-15  ID(N ,2 )  y-translation  boundary  condition  code 
Col.  16-20  ID(N,3)  z-translation  boundary  condition  code 

Col.  21-25  ID(N,4)  “-rotation  boundary  condition  code 

Col.  26-30  ID(N,5)  6-rotation  boundary  condition  code 

Col.  31-35  KN  Node  number  increment 

Note  that  an  unspecified  (ID  =  0)  degree  of  freedom  is  free  to  translate 
or  rotate  as  the  solution  dictates.  Delated  (ID  =  1)  degrees  of 
freedom  are  removed  from  the  final  set  of  equilibrium  equations. 

ID  =  -1  is  used  in  the  generation  of  boundary  condition  code  1.  Generation 
of  the  boundary  code  is  used  when  a  series  of  nodal  cards  all  have 
fixity  in  a  given  direction.  For  example,  a  flat  plate  lying  in  the 
x-y  plane  subjected  to  plane  stress  state  will  have  ID(N,3)  =  ID(N,4)  = 

ID(N,5)  =  1  for  all  the  nodes.  Rather  than  punching  "1"  in  column 

20,  25,  30  on  all  the  cards  it  is  possible  to  just  punch  "-1"  in  the 

column  19-20,  24-25,  29-30  of  the  first  nodal  card  and  enter  1  for  KN 

in  the  last  nodal  card.  The  program  will  set  ID(N,3)  =  ID(N,4)  =  ID(N,5)  =  -1 

on  all  of  the  intervening  cards.  A  code  of  -1  is  then  interpreted  in  the 

same  way  as  +1  (i.e.  fixed). 

Card  11  Material  Property  card  (8F10.0)  (one  for  each  1  ,pe  of  material) 

Col.  1-10  E( I )  El  longitudinal  Young’s  modulus 
Col.  11-20  PR( I )  VTL  Major  Poisson's  ratio 

Col.  21-30  E2(I)  Ex  transverse  Young's  modulus 

Col.  31-40  Gl( I )  GLT  shear  modulus  in  the  L-T  plane  of  the 
unidirectional  composite 
Col.  41-50  G2(I)  GT  shear  modulus 

Col.  51-60  VTT  minor  Poisson's  ratio 

Col.  61-70  WRANG(I)  ply  angle 

Note  that  a  different  ply  angle  is  considered  to  be  a  different  material. 

The  ply  angle  is  defined  by  the  angle  between  the  intrinsic  coordinates 
and  the  principal  material  coordinates. 

Card  12  Yield  strength  for  unidirectional  composite  (8FI0.0) 

(one  for  each  type  of  material) 

Col.  1-10  YLDX  tensile  yield  stress  in  longitudinal  direction 
Col.  11-20  YLDY  tensile  yield  stress  in  transverse  direction 
Col.  21-30  YLDS  shear  yield  stress  in  L-T  plane 

Col.  31-40  YLDXX  compressive  yield  stress  in  longitudinal  direction 
Col .  41-50  YLDYY  compressive  yield  stress  in  transverse  direction 

Omit  card  11  if  INCMAX  =  1. 


Card  13  Layer  information  card  (315)  (one  card  for  each  element) 

Col.  1-5  L  Element  number 
Col.  6-10  N LAYER  Number  of  layers 
Col.  11-15  KN  Generation  code 

If  KN  is  left  blank,  one  card  is  needed  for  each  element.  If  KN  is 
set  equal  to  1,  only  the  first  element  in  the  series  need  be  provided. 
The  other  will  be  set  equal  to  the  first  element. 

Card  14  Layer  property  set  (215,  Fio.4,  15)  (one  set  for  each  element) 


Col. 

1-5 

LN 

Layer  number 

Col. 

6-10 

MTYPE  Material  type  number 

Col. 

11-20 

THL 

Thickness  of  the  layer 

Col. 

21-25 

KN 

Generation  code 

KN  is  defined  the  same  way  as  in  Card  12.  If  KN  is  left  blank  one  set 
of  cards  is  needed  for  each  element.  If  KN  is  set  equal  to  one  at  the 
last  card  of  the  first  set,  then  only  the  first  set  is  needed  for  the 
first  element . 

Card  15  Cubic  spline  control  card  (15,  5X,  2F10.0) 

Col.  1-5  N  Number  of  segments  of  stress-strain  curve  to  be  fit 

Col.  6-10  Leave  blank 

Col.  11-20  E  The  initial  shear  modulus 

Col.  21-30  ES  The  last  shear  modulus 

Card  16  Discreet  values  from  the  shear  curve  (2F10.0)  (one  for 
each  station) 

Col.  1-10  F(I)  Value  of  shear  stress  at  station  I 

Col.  11-20  X( I )  The  corresponding  shear  strain 

Omit  Card  15  and  16  if  INCMAX  =  1. 

Card  17  Concentrated  load  card  (15,  5X,  3F01.4)  (one  for  each  load) 

Col.  1-5  ND  Node  number  where  load  applied 
Col.  6-10  Leave  blank 

Col.  11-20  IDIRN  Direction  of  the  applied  load 

1  for  x-direction 

2  for  y-direction 

3  for  z-direction 

Col.  21-20  FLOAD  Magnitude  of  the  applied  force,  positive  if  in 

the  positive  direction  of  the  axis,  negative  if  opposite 
to  the  direction  of  axis 

Card  18  Distributed  body  load  card  (15,  5X,  3F10.4,  15)  (one  for  each  element) 


Col. 

1-5 

L 

Element  number 

Col. 

6-10 

Leave  blank 

Col. 

11-20 

Px 

x-component  of 

the  body 

force 

per 

unit 

volume 

Col. 

21-30 

Py 

y-component  of 

the  body 

force 

per 

unit 

volume 

Col. 

31-40 

Pz 

z -component  of 

the  body 

force 

per 

unit 

volume 

Col. 

41-45 

KN 

generation  code 

Omit  card  18  if  LD  i  1. 
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Card  19  Edge  pressure  load  card  (215,  F10.4)  (one  for  each  load) 


Col. 

1-5 

L  Element  number 

Col. 

6-10 

IS IDE  Side  number 

Col. 

11-20 

PMLD  Edge  pressure  positive  if  in  the  same  direction  as 

the  outward  normal  of  the  edge  surface 

Omit 

Card 

19  if  LD  i  5. 

Card 

20 

Surface 

pressure  load  card  (15,  5X,  F10.4,  15)  (one  for  each  element 

Col. 

1-5 

L  Element  number 

Col. 

6-10 

Leave  blank 

Col. 

11-20 

PN  Surface  pressure,  positive  if  in  the  same  direction 

as  the  outward  normal  of  the  surface 

Col. 

21-25 

KN  Generation  code 

Card 

21 

Non-uniform  surface  load  set  (15,  1,  8F10.0)  (two  for  each  element) 

Col. 

1-5 

L  Element  number 

Card 

22 

(Continuation  of  Card  21) 

Col. 

1-10 

PU(L,  1) 

Col. 

11-20 

PU( L,  2) 

Col. 

21-30 

PU( L,  3) 

Col. 

31-40 

PU(L,  4)  Pressure  intensity  at  the  eight  nodes  of 

Col. 

41-50 

PU(L,  5)  the  element. 

Col. 

51-60 

PU(L,  6) 

Col. 

61-70 

PU(L,  7) 

Col. 

71-80 

PU( L,  8) 

Omit 

Card 

20  and 

21  if  LD  i  3. 

Note 

on  Region  Connectivity  Data  for  Mesh  Generation 

A  domain  is 
connected  to  one 

generally  modeled  using  several  quadrilateral  regions 
another  along  one  or  more  sides.  The  possibility  of  a 

common  boundary  between  two  regions  requires  that  certain  information 
be  provided.  The  determination  of  this  connectivity  data  is  best 
illustrated  through  an  example  as  the  four  region  body  in  Figure  59. 

The  £  n  coordinate  system  and  the  region  number  have  been  assigned.  The 
sides  of  each  region  are  indicated  by  the  number  1  to  4.  The  connectivity 
data  for  the  four  region  body  is  as  follows: 


Region 

1 

2 

3 

4 

1 

2 

3 

0 

0 

2 

4 

1 

0 

0 

3 

4 

0 

0 

1 

4 

2 

0 

0 

3 

The  first  line  of  data  states  that  side  one  of  region  one  is 
connected  to  region  two  and  that  side  two  of  region  one  is  connected 
to  region  three.  The  two  zero  values  indicate  that  sides  three  and 
four  of  region  one  are  not  connected  to  any  region.  There  is  one 
line  of  data  for  each  region. 


APPENDIX  B 


SRESET  LIST 

FILE  1 <TITLE="CH/TAPE".KINO=DISK#FILETYP£=?) 

CGMMQN/SET?/INC*NFLAG»YL0X(4)#YLQV(4).YLDSC4 >»  YLDXXC4 ) 
1#  YLQYY<4),SUMS< 100 * 8  *  3 > >K 0U NTC 1 00) ► 

2  Xf)UNTF(100)*NFAIL(  100.8  )*NYI  ELOUOO*8  ) 

3  *  SUM  SI  G<  190. 8. 3  )  ,  SUM  S  TN  (  1  00*  2  0  >.  SUM  EXC  1 00*  20> 
4.SUMSX2CI9C.20  ) 

CG-MON/Hl/ X(163)»Y(153)#Z(163  ). 10(163.5) 
CCNM0N/HS/E(4).PR(4).N0D(44.3) »MT)PE<0*44  >.TH(44  ). SE 
1(45.45  >.  PX<44  )#PY(44)*P2(  44).NLAYER(44)#THL 

2<t.44  )  .£2<4l*Gl<4)»  G2(4 ).PR 2( 4 ) .*R ANG< 4) 

C0HMQN/SPN/DDOC20O  I.NN3.E5 
C  I  HE  NS  I  ON  T  ITL  E  (1 3  ) 

3  111  EM  SION  F(5».S(320.80).LM(4C).FF(45).R(  320)  #S  1 1  9  )  .T  I 
1(9)  /  AA(2),FFF(6)#AX(136> 

CC*M0N/G/AMN(183.6  >.PM(44. 4) 

CGMHON/P/NSHELL*N9»NEOGEL»K01JNTy»KOUN 

CCMH3N/S0L/RR(302> 

CCMHON  /TEST/iNCMAX 
CCHHON/LO  AO/PNC  4  4  ).Ptj<  4  4*8  ) 

CCHM3N/L3A0GH/LD.NE3.NL0  AO.N£L*NELPL»NNP 
DATA  AA/-. 57735027. .57735027/ 

DATA  IPRC/IHF/ 

CATA  Vi  OR  01  /6HH  ONLY/ 

N9=e 
NFE  =8 
N  r '’  =  3  2  0 
hLC  =  c  -3 

RE AO (5. 901 )(T! TLE( I >.1=1.13) 

WFITE<6»90l)(TITLE( 1)*1=1.I3> 

901  FCRHAT(13A6) 

RE  At)  (  5  .  S  )  NNP.NEL.NMAT.NSHELL.NFLAG 
REA:H5.6)LD»NELPL.NE0GEL.NL0AD.  inch  ax.nelsp 
IFCNELPL.E3. 0)NELPL=NEL 
6  FCPHA1C6I5) 

IF(NFLAG.Ea.O)  NFL  AG  =  ? 

WFITE(5#20OP)MNP,NEl,NHAT*NEnGEL»M.0AD'NSHELl*NFLA  G 
1  .  INCH  A  X  .LO 

N  TOP  =5  *NN? 

1F(NNP.NE.0.AN0.NEL.NE-0)G0  TO  300 
C**  MESH  GENERATION 

CALL  MESH* (X.Y.7.N00.NNP.NEL.44.TITLE) 

IF(TITL'E(1  3).rQ.H0R01)  CALL  EXIT 
N I  OF  =  5  *NN? 

IF(NELPL.EQ.O)  NELPL=NEL 
GG  TO  305 
30C  CCNTINLE 

C*  *  I N?J  T  ( 3  D  4  L  C00P01NATFS 

call  Input  cx.y»7»nnp) 

C**  INR’JT  ELEMENT  CONNECTIVITY  OATA 
CC  521  H  =  t  .NEL 
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READ  C5,bl 5)N* (NQ0<N»1 )*I=l»3) 

ITE  (6>ol  5  )N»  (NOQ(N»  I  )»  1  =  1  #8  ) 

52C  CONTINUE 
3)5  CCNTINUE 
C*  *  INPUT  II)  COOES 

C  ALL  1 N ! 0(  10»NNP,NEQ) 

515  FC  PM  A  1 (9T5  > 

WnIT£(6»?OlC) 

WFITE(6f20Ot5)(N»(IO('l*l)*Isl*5)*X(N)*Y(N)*Z(N)*Nsl 
1»NNP  ) 

2 C C 5  F0PMAT(6I5,5X»3F13.3) 

9  32C  r  C  PM  A  T  ( 1  HI  ) 

4RI  T  E  (  6  >  20  2  )  ) 

WRITE (fi#2Q25)(N* CNOO(N»l  )  *  1=  l #8)*N=1»N£L  ) 

WRIT E (6, 2) 30  ) 

C**  INPUT  MATERIAL  DATA 

OE  5  I  =  1  *  N  M  4  T 

OtAO  (5»3>  E(I)*PR(I)#E2(I)*G1(I)»G2(I)*PR2(I)#WRANG 

1(1  ) 

1  Ft  INC  MAX.  NF.l  >R  E  AD ( 5»  3  )  YLQ  X(  I  )  *  YLO  Y(  I  )  *  Y  L  05C  I  >»YLOXX 
1(1) »  YL  0 YY(  1  ) 

hfil  TEC  5*20  35  )I  *  £(1)*PP(I)#E2(I  )»Gl(  I  )*  G  2( I  )*PR2(I) 

i»M«ANGCI>»  YL  OX  (  I  >#YLDY(1>*YLDS(I>'YLDXX(I) 

2,  YLDYYO  ) 

5  CONTINUE 
i  FCPM4T< «F1 Q. 0) 

4  FORMAT (1t»12E10.3) 

KN-0 

WRITE(6>204C) 

C**  INPUT  NUM3ER  OF  LAYER  FOR  EACH  ELEMENT 
SC  ir  M  =  1,NEL 

IFCKN.E  J.O  JREAD  (5 »23 >L» NLA  YE R(L >*KN 
If (KN.FD.O  ) GO  TO  27 

ML4YEF(M)=NLAYER(1  ) 

27  CONTINUE 

17  WFITE(6*204'5)M»NLAYER(M)fXN 
It-  FCMAT  (315  ) 

KN  =  ) 

Wf  IT  E(  6*20  50) 

CC  4?,j  M=1*NEL 
NPLY  =NL*YER (M) 

U  =  0. 

O*  IISPJT  MA  Tr  R  T  AL  TYPE  AND  THICKNESS  OF  EACH  LAYER 
OG  4 l C  L  4  =  1  >  NPL  Y 

1F(MN.E0.  ))REA0(5*12)LN#MTYPE(LN*M)^THL(LN»M)»KN 
;P(KN.E3. 0 >  GO  TO  411 
MlYPE(LA,M  )=MTYP£(LA»1) 

ThL(LA*M)=THL(LA*l) 

411  TL=TL«THL(LA*H > 

41  C  i*R  1  TE(b»?0  55)M#L4*MTYPE(LA*M)  *  THL (L  4*  M) *  KN 
TH(M  )  =  TL 
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42C  CONTINUE 
\  >  FC*MAT(2I3  »r10.4*15) 

IF(  I NC  MAX . E  C„  1  )  GO  TO  3030 
C*  *  CALL  CG3IC  SPlINE  FIT 

CALL  SPLINE (DOO  <1), 000(51 >  »DOO(  101>,0D0(15t)*NN3»E5> 
3030  CONTINUE 
C**  INPUT  LG  AO  TNG 

CALL  I NLOAD(R*PX,PY  »?2  ,PM,PN,FU»NE0GEL> 

M  t  A  N  0  =  0 
LlM=NPE*5 
KCUNT Y=0 
CC  ?05  M  =1 ,  NEL 
KC'JN  TF  CM  )  =  3 
(COUNT  (M)=0 
JO  ?  )5  L  =1  ,  NPL  Y 
SLMSTNC  3*L  )=0. 

SlMtX(M#L)=C. 

3l*S<2  (*«»L  )  =0. 

OC  005  1=1,3 
3li^3lG(M,L,  I  )=  0  . 
nS  SLM5C  M,L»  I  )  =0. 

GG  60C  INC=1,1NCMAX 
CC  605  1=4,6 
505  FFFC I)=0. 

IFUNC.GE.2. AND.KCJUNTY.NE.I  >GC  TQ  505 
CC  U  11=1, NEQ 
RFC  11  )=H(I 1  ) 

ca  li  ji=i,noc 

11  S(Il,Jl)  =  0. 

CC  5  OC  M=1  , NEL 
CC  5)1  I V  =  1  , NP E 
I i=5*(IV-l ) 

N 1 =N0 C ( M, 1 V  ) 

CC  301  JY=1,5 

ij=ii*jv 

501  LN(IJ)=I0(N1»JV  ) 

~i  FCKNA  T(  201  5) 

C  * «  CALCULATE  5ANDWI0TH  AN  0  ELEMENT  STIFFNESS 

:all  eancal  (toand,l;i,.npf> 

CALL  ELEMN I (M,FF  ) 

CC  A  1 5  LL=1,LIM 
r=LM(LL) 

IFU.LE.O)  GO  TO  A15 
3f(  I  )  =  Fr  (  I  HFF(LL) 

CC  A0-;  MM  =  1  , L  1  M 
J  SLN  C  HM  )"I *1 
IF(J.LE.O)  GO  TO  AOO 
S(T,J)  =  S<I,J>*SECLL,MM) 

400  CLN TINGE 
A 1 5  CGN7INUC 
50.  CONTINUE 


ta  FC.SMATUX,5H  NEQ=,  15,  5X,  7H  M3AN0=*I5) 

**  DISPLACEMENT  CACCLATION 

CALL  3  ANSGLt i,M9AND,NE9,Rfi  >S*NDR,N0C) 

CALL  3AN'$JL(2»M'TANO,NEQ,RP,S,NDR,N0C) 

A flTE(&,?CcC) 

505  CCNTINLE 

CC  7  r /  N -1 , NNP 
DC  7  7  6  j= L , c 
F  t  J ) =0  . 

1 i=IOtN, J) 

IF(II.EO.O)  GO  ra  77 5 
Ft J) 1 ) 

776  CONTINUE 

CC  775  1^,3 
7  73  FFF< I  )  =F  t  I  ) 

FFF  ( <•  )=AMNtN,l)*Ft4)4AMNtN»4>*Ft5) 

FFF<5 >=AMNtN,2)*F(4)* A VN C N , 5 ) *F ( 5  ) 
FfF(o)=AMNCN,3)*F(4)4AMN(N»6)*FC5) 
IFdNC.EO.  l)Wt?ITE(b.0  02)N,(FFF(J6),Ib=l,6) 

777  CONTINUE 

IF(INCHAX.“C.ir»IHlTE(6»2070) 

I  Ft  INC.  EJ.  2  )WF  I  T  C  t  6,  2  OPO  ) 

**  S  1  f'E  3  S  CAC  UL  A  I  ION 
OG  30C  M  =  1  , NEL 


CALL  STRESS  t M ,N  ELSP  > 


510  CCUTINUE 

DC  1777  N  =  1,W.NP 
CC  1  776  J= 1  *  5 

F  t J) =C  . 

I  1-IKN,  J) 

IFtlI.El.O)  <10  TO  1776 
Ft  DIRECT!  ) 

177=  CCN  T  i  nj OF 

<  tN  )  =  X  t  N  )4  F  f  1  ) 

Y  tN  )  =r  <N  >t  F(  ?> 

Zt.N  )=i  tN)4  r(  3) 

1  77  7  CONTINUE 

90  2  FCP*  A  HI  X,  I  5,  3T  U.4,4X,  3F  14.4  ) 


SCO  CONTINUE 

ICC  FCPMATtlHl  1  X,  37  HC  ON  T  ft  0  L  I  N  F  0  R  M  A  T  I 

1  C  N  7/  1X»21HNUM3EP  OF  NODAL  POINT  29 

2tlH.l  1 H  =  15  /  IX, ldHNUM  3C  F  OF  ELEN  tN f  3? 

JtH.)  1H=  !>  /  l  X*  1  6  HN'J M  G£  F  OF 

4  MTEPIAL  32tlH.)  1H=  15/  IX/32HNU. 

5  OF  ELEMENT  WITH  EDGE  LQAOING  i«CH.>lH  =  15  /  U 

6,24HNG.  OF  CONCENTRATED  LUAO  26tlH.  )  IH-  IS  / 

7  1  X  ,6  HN  5  H  ELL  44(1H.  )  lH  = 

e  Ii  /  1 X, 2  7  H  Eti.,  (  FOR  SHELL  ANALYSIS 

?  /  ix, pan  :a.,  1  fo«  plate  analysis 


ix, lhnflag  4 5 (  1 H • ) 

/  IX,??H  F5.»  1  HILL  C^ITEPIA 
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1X.25H  FQ  •  *  l  TSAI  "MU 

/  1  X  ,  1  P HH  AX. 

1 H  =  i5  / 


4  CRITERIA  /  IX, 1 

5  INCREMENT  NO-  32(iH.)  1H=  IS 

O.9HL0AD  fY°E  41(1H.)  IH=  I*  /  IX  44 

?  cistr nurro  loao  in  n  z  oirecthn  / 

6  £0..  2  NORMAL  LOAO  /  IX  36H  £0.*  3 

)  £:STh!«UTEO  LO  AC  /  IX  ?6H  EU..  4  CO 

1  LOAO  /  IX  16 H  E9.,  5  E3GE  LOAO  / 

2  02  5  FCfiMA  1(1 3* 5 x# 61 5  ) 

200  FCi'lAK  //  IX  (IMMATERIAL  10X  IHE  IX  2H 

1  dX  2 1*31  6  X  »  2HG2  7 X  3HPR?  5X  5HWRANG 

2  £  X  4HYL0Y  6X  4HYL0S  5X  5HYL0XX  5X  5HYL0YY 

3  ?n  'IUMOro  ) 

2;33  rC4MATCl5»rX»l?ElJ.3) 


is  /  IX 
4  4H  2  9.*  1 
IX  20H 

3  NONUNIFORM 
C  ONCEN  Tk  A  TED 
) 


2  HE  2 


8X  2HPR 
6  X  4HY13X 
/  IX 


2C2C  rC’^H 
1 
2 

3  2  3 

2  C  1C  FORMAT ( 
1 


1 M I  3  4  H  ...  2  NOOC  THIC<  SHELL  ELEMENT  OATA 

//  IX  ?HELEMENT  isx  I2HC0NNECTWI  TY 

/  1 X*  5h  NQ.  48h  l 

4  5  6  7  r-  ) 

////  L  X  ,2  GHGENER A  TEO  SODA  L  OATA 

////  1  X,  16HE  9  UA  TIO  N  NUM3E9S 

////  1  X  »  3 1  HNODE  DEGREE  OF  F  0  £  E  0  0  M 


1  X.2  3  HNOOAL  POINT  CCOROlNAfES 


4  /  IX  3  4HNUM3FR  X 

5  9  X  IRX12X  1HY  1 2  X  1H7  ) 

2C4C  FCRMATC////  IX  22  HELEMEN  T 

1  ifiH  NO.  NO-  /*> 

2  3  4  3  FOHAT  (IX,  15,5 X,  1 3.7X,  12) 

2C5C  FCRM  A  T ( /// /  IX  39  HELEMENT 

l  THICKNESS  /  IX  16H  NO. 


LAYER 


layer 

NQ. 


ALPHA  OATA 


MATERIAL 

/,  ) 


2  05  5  FORMA T( IX. 15,5 X, I3*7X, I  3, a X, F 10. 4, 1  5) 


’05C  FORMAT  UH1 

1  T  S  /  P 

2  // 

3  IX  4CH  7  - 


IX  40 H 
IX  40H  0 
IX  4)M 


MOO 
T  A  T 
N09c 

X- 


i  40HNLM3ER  TRANSLATION 
.  40HNSLAT! ON  ROTATION 

ROTATION 

FORM  A  1 ( IH1 1 X  45H...S-N30' 


20  7  C 


1  COM  °  0  TED  . 

2  EHSK-XX 

3  SIG-XY 


IX  4  5  H  .  •  •  3-NOOE  THICK 
A  T l X  PHCENTROIO  //  IX 
“X  IMS  49HIG-YY 


Y  IX  40H 

/  lx 

TRANSLATION  TR A  IX 

ROTATION  IX  40H 

/  ) 

SHELL  ELEMENT  STRESS 
l 4H EL  tMENT  LAYER  IPX 
S IG'22 


5  21HIG-Y2 

7  ) 

2  C  P  0  FOt*M  4  TMH1  IX  43H  INC  GL.  NO.  LAYER  YlELO 

1  SIGX  SIGT51H  SIGXY  SI 

2  S2  S1216H  STRAIN/) 

CALL  exit 

FRO 

SUBROUTINE  ELEMNT(MM.F) 

CCMM0N/SET2/IN'C,NFLAG»YL0X(4)  »YLDY(4),YLDS(4),''L0XX(4) 
1.  VLO  YY(4  )^SUM  S(100»  E»  3)*K0’JNT(  IOC)  , 


SliWX 


YlELO 

SI 


•  • 
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2  hOUNTc( 130  l.NFAILCtOO.S  ). N Y 1 ELD( 1 0 0.  8) 

3#  SMb  I  1(100  *3#  3  ).SUmSTN<1  ?0.2  0>.SUNEX(10'),'*Ci> 
CCMMON/YIEL0/LL.MT»EL.ET*GLT*¥TL»YTT.G23» G1 3»  ALl* 1PC 
CCMM0N/Hl/X(16  3).Y(i&3).Z(lS3l»I0C163#5) 

COMMON  /H?/XL(9  )  .  YL  <  9  >  »  ZL  (  9  )  It  9>»tf 2(  3>.V  3<  9)  * 

1  J4C(  3.  3)  .N(9>.NS(9j.NT(9)  .01  <9  >»0  2(9  >.  0  3 

<?(9).W1(9).  W?(9>.W3(9) 

CCMM0N/HP2/TT.U11.U21. U3t. U1 2 . U22, U 32  .  I 
CCMM0N/H3/E(4).Pft(4).N00(4<.»9  >.MTYPE( 3.44).TH(44).$E 
U  45.45).  PX(44  ) »  P  Y (4  4  ).  Pi  (44  )»NLAYER(4<*).  THL 

2  (  8  #4  4  )  .E2(4).G1(4),G2(4)>PK2(4).WRANG(4) 

3  .Q(5.6).H(9> 

COMMON/ SPN/OD0(2OO  >»NN3.E5 
CCMMON/G/AMNd  6  3.  6  )*PM  (  4  4.4} 

C  1ME  N;  I  ON  0(6.45  ).33(o.45>.KKK(6).E3CS.45).D3(6.*5).EE 
Ufc.6>.M(2>  .FX(45)»FY(45)*FZ(45)#F(45)»AA(2)*DE(b.6) 

2.CD( 6.6).0EB(6.45)  *0E3l(6.45>  »FF2(5> 

3  .P0(9).FXX(9  }.  FYY(9  ) *  F22(  9)  .E£SAV(6.6> 

4»CESAtl(6»6 )»0QSAV(6»6) 
CGMrf3N/P/NSHELL.N9.NEQG£L»KCUNTY.KaUU 
CCMMQN/L0a0/PN(44  ).HU(44,o> 

C  CMMON/LQA  OGH/LO 
COMMON/ K/K l  . \2*  K3.K4.K5 
CATA  H/1..1./ 

SATA  AA/"*«  5223502  /  *  •‘5  77  35022/ 

SEAL  N.NS.  NT.JAC.Ntil2.Nw2l 

NNN=INC-1 

N45  =  N9*5 

30  100  11=1.3 

XL(!1 >=X<NOO<MM.ll>> 

YLC 1 1 )=  Y(M  UD  (MM . 1 1  )  ) 

2  L (  1 1  )=Z(NOO(HM*  II  )) 

PC<  I  l  )  =P1J(  MM. II  ) 

ICC  CONTINUE 

105  FfRMAIOFl  0.0) 

>L(9)s-(XL(1  )«XLC3)«XLC5  >•*  XLC  7  >  )/4.  *(  XL(  2  >♦  XL  (  4  >♦  XL  <  6  > 
l*XL(  3)  )/?. 

YU  >>*-<  YL  (t  >♦  YH3  WLC5  WH  7))/4.*(  YLC  2)»YL<  *>*YLC6) 
lML<3»/2. 

2L<  9)  =  -(2L(  1  )4ZL(3)-»ZL<5  )-»ZL(  7  )  >/4.a(2L(2  )*ZL<4  )«ZL(6  > 
1 ♦ ZL( 8 )  )/2. 

PC (9  M-(P3(1  )4P0(3>4P0(5)*P0(  7  )  )/*.♦( PQ<  2>*P3(  A  >4PJ<6  ) 

1  ♦  F 0  (  3  )  )/2. 

GC  TO  ( 1 .2  »  6 .6  »  6  )»  LO 
1  CONTINUE 
PTX=PX(MM) 

P1Y=P Y(MM> 

PTZ=?Z(MM  ) 

GO  TO  6 
?  continue 

F*L 1AC  =  ?HCMM  ) 
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e.  continue 

CC  2s  «  =  i»AI45 
F(K) =0. 

F  KK  )=0. 

FUIO-O. 

FiCK>  =  0* 

DC  26  L  =  l*N45 
26  SF(K#L>=0. 

CC  2  5  X=!  »  b 
0C  25  L  =  1  *  6 
J£CK,L  >=P. 

CC(K,L  )=  0. 

?5  EE<K*L)=0. 

NPL  Y  =!\iLA  YE  ><  (  MM  ) 

T  T  =  T  Ht  rtrt ) 

H(l  )  -■*  TT/2  • 

0Q  35C  I  =  1 »  NPL  Y 
35C  H(l4i  )=  H(  l  )4THL( 

00  400  LL  =  1  • MPL  V 
NFA1L(MM#LH  =  0 
\Yl£LC(  4H*Ll>=0 
MI=iTYPr:(LL»»1M  > 

1F(IMC.£<3.  1  )  GO  TO  6001 
SU2  =  SUHSTN(M*,LL) 

Xl=A9S<SNJ  ?) 

CALL  FUNC  T ( 009 ( l J*  000(5  l  >*090  (t'H  ),  000(151  ),Yl,NN3,Xl 
1*YP,  Y  F  P  »  f  5  ) 

GH  HT  >=YP 
u  E  {  i  T  )  •=  YP 
6  CO  1  CCNTIMJE 

Ul»i«MNij(MT  )*  3 .141 5*2654/1  GO. 

TU*(H<LL*  1  Y“H<LL))/TT 

TI2=(  (  2.<M(LL*  l  >/TT>**2-<  2-*HtLL  >/TT  )**2  >/4. 

TI J=C(2.*H(LL4l )/TT)**3-(2.*H(LL)/TT»**3) /6. 

£L=E< M  T) 

£  T  =  £  2CMT ) 

V  IL  =  PF (MT ) 

4  U-PR2(nr  ) 

GLT=ji(MT  ) 

G  T  =  G2 (MT  ) 

IFCET.E3.!).  >ET=EL 
IFtGU.EQ.0.)GLT  =  EL/(2.*(  1.4YTL  )) 

I F  (  V  Y  T  .  £0.  0  .  )  VTT  =  VTl 
IF<GT.€3.. i.  )GT  =  ET/C  ?.*U  ,*V7T)) 

GE  3  =  GlT 
G 1 3  -  G  T 

IFd'JC.EQ.tYGU  TO  3)2 
IPC  =  0 

CALL  YIE»-9(MM) 

3CE  CCNTINUE 

CALL  JMAT ( 9  ) 


zeci. j)=zeci.j>4  qci»j>*tti 
0  £  (  I  »J)=DECI»J)tQCI*J)*rT2 
OOCI#  J)=C>0  Cl  »J)*QCI  *J)*TTS 
52C  CCNTIKUE 

IFCNFAlL(MM»LL).eQ.l.ANy.K‘]UMTFCHM).Ea.2)  CALL  FXIT 
ACC  CONTINUE 

CC  5 3C  1=1  .f 
C  u  5  ?  vj  J  -1  #  S 
I£i4VC!»J)=£E(l.J) 

C E S A V< I* J)=0EC I#  J) 

:CSA  f (!* J) =  00  ( 1  *  J) 

5  *  C  CCNTTMJE 

IFC NShfLL. M£ - 0 >  GO  TJ  54 
CALL  SUPVGC(TT) 
jC  TO  59 
54  CCNT !NL€ 

OC  5 o  !=l»Ni 
VIC  i  )  =  c . 

V2C 1 ) =0. 

5c  V  3 C I ) =  T T 
59  CC.NTIMJE 
CC  50  I=t.S 
OC  51  J=1.N45 
3SC  I  *  u  )=  0. 

SC  3(1,J)=0. 

CC  200  KK  =  1 .2 
JG  20C  11=1.2 
CC  20C  JJ=1#2 
SS=A4(  II  ) 

T 1  =  A  A  (  u  J  ) 

2H  =  AACK'C) 

*1  =  «<C  II  ) 

W„=rfC JJ  ) 

HK=W(KK> 

CALL  SHAPE ( 5S. T1 » ZK.OETJAC  ) 

OG  5  A  C  1=  l  .  f. 

CC  540  J=!  .5 

EEC  I  . J ) =EE SAVC1 . J) 

Ot(l»J)=OcSAVCI  »J) 

CCC I. J )=3DS*VC I. J) 

540  CONTINUE 

I  F  CSS  HELL.  Z  3.0)  CALL  ZTMNCEE.DE.OD) 

X  MIL  l  =  H  I*WJ*WK*QET  JAC 
?NN=0. 

V  V 1 X  =  C  • 

VV2Y=0  . 

W3Z  =  C  . 

CC  1  S?  Trl  ,N  ? 

VV 1=  VIC  I  )  / T  T 


=  )/TT 
W  5=7  3(1  J/TT 
V»1X  = vVIX* VVl*N( 1 ) 

VY?Y  =  \lV2Y«  V  V2*  N  (  I) 

Y»5 2-1V17*  VV3*N(I) 

IF  (LO.tO.  3  )  °NN=PNN*P3(1 )*N(I ) 
15  2  CCNTINUE 

JG  ■>»  I  =  1 »  N 9 

I  F  ( N  5r*ELL<>  Nr  •  0  )  GO  TO  61 
U  1 1 =  01(1 ) 

L21=02( J  ) 

1)21  =  02(1  > 
tJl?=-Wl(  I) 
l22  =  -ii2(  I ) 

'J  2  2=  “*»  3(1  ) 

GC  rn  57 

61  CONTINUE 
U  11  =  0  . 

621=1  . 

■J  2 1  =  u  • 

Gi2  =  1. 

U22 -  0 • 

G3;2=o. 

57  CONTINUE 

ifd.Ea.9)  r-n  tc  55 
*=N03  (  (N,  1  ) 

*HN(K,1  )=Ull 
AN(K*2)=J  21 
APN(X*3  >=U31 
AMN(K*4)=-U12 
*PN  (  'X  »  5  )=  “‘J  ?? 

AM(K»o  >=-U32 

55  CONTINUE 

C«Ll  8MAT(3/30) 

GC  T  0  (  i  6  » 1  2  *  1  3 »  16#l6)»Lr) 

12  CONTINUE 
PTX=?Ni.0A3*VVlX 
PT  Y  =  0.NLQA3*  W2Y 
F1/  =  PNL0A0*  V  VS7 
GG  IJ  16 

13  ?TX*PNN*mX 
P  1Y  =  P.nN*'7Y?Y 
P!2=°Ni'WV3Z 

16  CCMTiNUE 

FXK1)=N(I  ) *®T  X 

FY(K?)=N( I  )*PTY 

F2CX3  )  =  N ( I  )*P T2 

r>(K4)=NU  )*U11  *TT*HK*PTX*.  5 

FX(X5  )  =  N  (  T  )*Ul?*TT*W*(*PTX*-5 

f  Y(X4  )  =  N(1  )*U21*  TT*.(*(*pTY*.5 

F  I  (  K3  )  =  N( 1  ) *U2P* TT*WK*PT Y*. 5 


Fi(K4)=MtI  )*,J3!*TT*rtK*f>TZ*.5 
FI  i  *  5  )  =WM  )  *  U3  2*TT*  WK  *?  TZ*.5 

5  ?  CCNTJMJE. 

JO  60  5=1,6 
OQ  60  L=1»N45 
IE(K»l  )  =  ':. 

3E3 ( K*L  )  =3  • 

3EU(K»L)5y. 

Oe( 5  »l  )- 0. 

CC  b  0  M-t  »  o 

E  E(5  ,1  )=G8  <K»L  MEE<K,M  >*8CM,L) 
3EB(K»L)=0£«(K»L)40£(K»*0*3a<  ) 

CE31  (K*L)  =  0r91{  K»LMOE(K#M)*gtM*L> 

6'-  C  S  (K  »  L  )=  3  3  (  K.» L  ) ♦  DD  (  K*  M  >..* 39  CM*  L  ) 

.0  0  7  0  K=l»N45 

F< K )=F  U  >«  (FX(K  MFY(K)4F?(X))*XMUL1 
:c  n  l=i,na5 

OIM=0. 
i3l*l  =  C. 

CUM^=  j  . 

3  0*3=0. 

00  65  M  =  1  *  6 

CLfU  =  tUMl«  9  <9#K  )*E3«.M»  L> 

JlM2=JLMZ*9M,K>*t)£8(M,L) 

OOM  5  =  0UMT*  a3<M#K  )*0E-3l<M.L) 

6  5  CL't=  JU.i  +  rH  <M.  L) 

1  C  5t(  <»L  )  =SE  <  K,L  >♦  <00.9*  36"  l  ♦O'JM  2«  O'JM  3)*XMoL  1 
ro  f  ala  =  tomla*  xmul  i 

20  C  CONTINUE 

CC  512  I=l,N4  5 
FC1  )  —  F  C  I  >/TT 
012  CENT  I N  U  E 

TOTAL <=TOr ALA/TT 

IFlMECOEL.Nt.J)  CALL  L OCAL (F XX p F V Y ,FZ 7 *A A , W 

1 ,  fH  > 

CC  509  1=1 *5 
5'!  9  FF  ■?{  I)=0. 

JG  51 5  I  =  l »  M  9 
F  F ?  (  1  )  =F  XX  (  ?  ) 

FF2C2 )=FYX  (I  ) 

FF2(  3)  =  FZ7 (  I  ) 

OC  5  19  J  =  1  > 5 
I  A  =  (  1-1  )*•»♦  J 
515  FC!4  )=F(I4  MFF?<  jJ 
3ETUPN 
END 

SL860LTINE  SHAP  £  { S ,  T,  2K,  Of  TJA  C  ) 

CCMH  :)N/H?/XL(‘1>,  YL<9)  #7LC  9)  Zc  9>*  vs<?  >» 

1  JAC(3*3),N(9 )*NS(9)»VT(9) 

DIMENSION  JACI ( 3,3) 

CH.-MSION  Sl(5),TIt8) 


cm  si/-i.><j.*i.*i.#i.*o.*-i.»-i./ 
cat  a  ri/-;.*-i.*-i.*c.*t.*i.*i.*c./ 

*£4*_  jAU 

SEAL  N*N5*NT*JAC 

)C  11  1=1*9 

CliM  1  =1  S I  C  I  )*  s 

31^2=1. ♦  ti  c  I  >*  r 

If  (  1  .tJ.’.O*.!  .EQ.b)  GQ  TQ  11 

IF  (  t.l'j.4.  js.i.ca.ft)  G'l  TO  l7 

)LN3  =  S * 31  (  I )  *T*T  KD-l. 

N  ( 1  )  =  C'JMlOUM2*DUM3/4. 

IS  SC  I  >  =  (O'JN  1  *i3(JM2aOUH2*OUH3)*3  IC  !  )/<►. 

M(l)  =  (3UM1  *0|JH?*0UNl  *0  JM  3)*TICi  )M. 

GC  7.3  11 
13  :tN4=  1  ,-i* *2 

M(I  ) =CoM2*  00K4 /2. 

NSd  )=-s*o on? 

MCI >  —  T  I  ( I  ) *DUM  4/ 2  • 

GC  TO  11 
17  CLM5  =  1.-T**? 

ISCI  )=lUH1*DUM5/' 2. 

NSC  I  )  =  0JM5*SIC  I  )/2. 

NlCl  )=-T*3UMl 
11  CONTINUE 

N  m  =< i.-$ **?)*< l.-T**  2) 

NSC  7)=-2.*S*Cl.-T**2) 

MC9  )  =  "2  •*  T  *C1. -$**?» 

DC  15  1=1*3 
CO  15  J  =  l*  3 
lb  J«C(I*J)=0. 

3C  20  I  =1*9 

JACtl*l>=JAC(l*l>-*NS<l>*XL<!MNS(I>*ZK*Vl<I>/2. 

JAC<l*2)=JAC(l*2)-MS<I)*YLCIM*SCl)*Z\*</2(I>/2. 

JACC  1*3  )=JACC1*  3  )  ■♦  N  S  C  I  )  *  ZLC  I  )  AN  SC  I)  *  ZK*  V  3C  I  )  1 2  . 

J  AC  C  2  *  1  )  =J  AC  (2  *1  )*  NTC  1  )*XLCI)*NTCI  >  *Z  K*V  1CI  )/?• 

J  ACC?*  2)=  J  ACC2*2)4NTC  I  )*YcCI)«NTCI  )*ZK*</2Cl  )/2. 
JACCZ*3)  =  JACC?*3>«NTCI)*ZL<T)«NTCJ>*ZX*v3CI>/2* 

JACC  5  *l>  J ACC  3*1  )«NCI )*V1CI  K2. 

JACC  3*  2)=J ACC  3. 2 )«NCI  >**  2( I )/2. 

JACC  3*3>  =  JACC3*3)«NCI)*V3CI>/2. 

20  CCNTINUE 

0  £  T  JAC  =JAC  C  1*1  ) * JACC 2*2 >*J ACC  3*  3)  ♦ j AC ( ?*  1 ) * JACC 3 

1  *  E ) * J  AC  C 1  »  3)  ♦ JACC 3 *1)*J ACC  1*2)* JACC 2* 3  I 

2“ JACC 1*  3 )  *  J  AC  (  3*1)*JACC2»?)  -JAC(i*2)*JAC(2*l  > 

3*  JAC  C  3*  3  >  -JAC(2*3)*JAC(3*2)*JACCl*l) 

JACK  1*1)=  (  JACC  2*  2)*  JACC  3  *3) -JACC  2  *3  )*JAC(  3*  >)  ) 

1/CET JAC 

JAC1 C2*1)  =  “CJACC  2*1)*JAC( 3*  3)-JACC?*3)*JACC  3*1)) 
l/DET JAC 

jACIC  2*1  )=  C  JACC  2*  1)*  JACC  3*2) -JACC  2*  2)*  JACC  3*  1  )  ) 
l  /  Cf  T  JAC 


J£CI< 1#2)=-<JAC(1»?)*JAC<3* 3)“JAC<1 »3)*JAC( 3*2)) 
1/C'ITJAC 

JACI<2»?>=  (  j&C  (.  1#1)*JAC(3#3)-JAC(1  #3  >*JAC(3#1  >  > 

i/crr j«c 

JACI < 3#2)=-< JAC(l#l>*JAC<  3#  21-JACC 1#2>*JAC( 3#1> ) 
i/atrjAC 

J  AC  I < l #  5  )  =  < JAC< l» ?)* JAC (2#  31- JACCl * i ) * JAC<  2#  2) ) 
l/o:  T  JAC 

jACI (  2  #  3)=“(J4C(  l#l)*JAC(2#3)-JAC(l »  3)*JAC(2#  1) ) 

1/CtT JAC 

JACI ( 3*  5)=  < JAC( l»l  )*JAC( 2# 2) -JAC Cl #?>*JAC( ?#1> ) 
1/ufTJAC 

cc  n  i=i # 3 

CC  5;  J  =1  *  3 
3  C  JAC(I»  J)=JACKI  »  J) 
nETU^N 
1  NO 

SLOHOUT  INt  SURF  AC  (XL#  YL>  ZL #Yl #V 2#v  5 #THI  CKL) 

Cl  HZ  MS  ION  XKl  )»  YL<1>#7L(1  )#  Y1(  l  )#Y2<  1>#V  3<  1)  #KPTS(9 
1#  A) 

DATA  KPTS/  2  #  3  #  A  #  5#  6  »  9  #8  #9  »  4  #  1#2#3#4#5 

1  #  e  #  ?  »  8  r  ?  *  8#  9  # 2  #9  #4  #  5# 6  #  7  # 6# 

?#2#3#4#5*6*7»3»9/ 

OQ  l  I =L # 9 
KF1  =KP  T5(  I  #  1  ) 

KF?=«PTS( 1 #2  ) 

*r  3=-(P  TS(I  #  *) 

KP4»SPTS<1 #4) 

A  L l  =  XL ( K3t  )-XL(KP2  ) 

API  =YL(KrM  J-YLCKP2) 

AN!=ZL(KPl )-7L(KP2> 

AL2-XL(I(P3  )-XLC(04) 

A *2= Yeti's )”Yl(NP4) 

AN2  =  ZLaKP i 1” ZLCNP4 ) 

1L  3=AM*AM2-AM?*AN1 
A«3=AL2*ANl-ALi*AN2 
AN3=ALI*AN2-AL2*AM1 
7=39RT(AL3**2«4M3**7«AN3**?> 

V\(1=AL3/Z 
¥<*2  =  AH3 n 
V»3=AN  Hi 

V  V 1  =  ’/Vt*THlCKL 
Yv2=¥v7*THlCXL 

3=¥rf 3*ThICKl 

¥1(1  >  2  V  Yt 
¥2(1 >  =  V¥’ 

V  2 C  I  )=  YV  3 
1  CCNUNUC 

RETURN 

£NO 

SL39GLTIN"  0  AN 3GL ( KKK # 1 9 ANO #N C Q 'H » A K# NOR # NOC > 


DIMENSION!  R(1 )*  AKCNOR*NOC> 

NFS=Nitf-’ 
if*  =N"3 

1F(«KK.EQ.2)  iO  TO  20  0 

CC  12;  N  =  1*NRS 

M=N"l 

M  F  =  M IN  OC 18  ANO»N* “M  ) 

P1VUT  =  AK('!*1  ) 

0  0  t?C  L  =2  »  M  ft 
CF=A*(N*l> /PItfUT 
1=M4L 
J  =0 

OC  lie  K=L  *  MR 

J  =  J4  i 

UC  AK(I»J}=AK(1*J)-CP*AKCN»K) 

12  C  A  K  (  M  »  L  )=  CP 
GC  n  4 00 

20G  OC  22 C  N=1 »NRS 
M=N*1 

MR=MIN0(I8  ANO»NR-M  ) 

C  F  =  F  <  N  ) 

R(N)=CP/AKCN#i  ) 

OG  22C  L=?*MR 
1  =  M  ♦  L 

2 20  R  Cl  )  =MI  )-4*<N.L)*CP 
iUNH  >  =  FC’i7  )/AMNR*l) 

OC  320  1  =  1  # NR5 

-I 

K  -  N  •  l 

MF  =  ^IN0(I8AN0»NF"M  ) 

OG  32C  K=?*MR 
L  =  M-»  X 

320  R  CN  )  =R  (N>- AXC.M*  K  )*RCL  ) 

4 0 C  RETURN 
END 

3L8R00mC  8ANCAL(H8AN0»LM»NP£  ) 

JlMF.^SION  LHC  l ) 

H 1m=1 000)0 

M  A  X  =  0 

L1,«  =  NPE*5 

OC  IOC  L  =1  *L  IH 

IFCHCL  >.£Q.O)  00  TO  100 

IFUH(L).GT.MS’OMAX=L*<l  > 

IFCLH(L).LT.H1N)MIN=LHCL) 

ccnttnue 

NO  IF  =K AX-MI NA1 

IFCN TlF.uT ./ 3AN0  )  M3AND=N0IF 

RETURN 

END 

SLIPOUTIn:  STRE ibCM-WNE LSP) 

CC'1MJN/i"T2  /  INC*  NrLA0*YL0X(4)  »  f  LU  Y  (  H)  >U0${  4)  ,YLUXa(4) 
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1*  YL0YY(4)*5'JmS(1M*  fc*  3>*K0'JNT(  100)* 

2  KilU3TFCl03)»NFAlL<l3n»?  WNYIELOCIOO**  > 

5»SU**1GC19Q  *-*  i  >*  SUMS  Tim  <  100*20*  3UM£X(  100*  20) 

4  *SiHSX2n>>'>»2‘' ) 

CG*»M  )N/Y1EL0/LL#MT*EL»ET*GLT*  VTL*VTT*G23*  G1  3*  AL1 *  IPC 
CCM10N/H1/XC  163  >*Y(  163).  Z(1&*)*I  0(163*5) 

CCMM  JN/H’/ XLC9  )*  YL(9)*ZL(? )*VlC9)*tf2( 9)*tf  3C9)  * 
i  JAC( 3*3) *H(9 )*MS<9>.MT(9)  *0  l(  9)  *02(  9)*  03 

?(5<*41C9>*  K?(9)*W3(<5) 

CCMM0\/H2_>/tt*UI  1*  U 21  * U  3 1  * U  1?  *  U22  *‘J 3?  *1 
COM3N/H3/ E(4)»PR(4)> MOO  (  44*3)  #MTYPE< 3*44  )*fH(44  )*SE 
1(  45  *4  5  >*  PX(44)*PY(44)*PZC44)*NLAYEF(44)*THL 

2C  8*4  4  )  *  £2  <4  >  *  G  1  <  4  )*  G2  (  4  )*  PP  2 ( 4  )*WRAmG(4) 

3  *Q(i*6)*H(9) 

COMMON/ SPN /OOOC  2  00)*NN 3*  E5 

C  CMMON/G/A  MNClo  3*6  )* PMC 4 4, 4) 

COMMON / T“ 5  T /1NCMAX 
CCMM ON/iC/Xl  *K2*  K  3*K4*K3 

C1M£MS13M  3(6*45  )*  33  (  fa  ,  45  )  *  KKKC  b  )*  E3C  fc*  4  5  )  *Dti  C  b  *45  )  *E  C 
1(«*6)*N(2)  *0(45)  *AA  AC  1  )*  SIGN(6) 

?*SIGM  (6>*STN(5  >*  STLC6)*  ST(6)»S!C,(6  )*STR(  3  )* 

3  S I GG ( 3 )*X I (9)*  YIC9 )*  2  I < 9)  *SI ( ?) *TI ( 9)  *STL2 

4C6)*31N2(o)*SIG2(6) 

COMMON /P/0  SHELL  *N9*NE0GEL*K0UMY*K0UN 
CCMM JN/SGL/UC  302  ) 

DATA  S1/-1 .*0.*1 •*1.*I.*0. *-l.*-l. *C./ 

DATA  11/-1  . *-!.*-!. *0.*1.*1.* 1**0. *0./ 

OATA  W/l.*l./ 

REAL  N*MS*NT*J*C*NU12*'JU21 
N  45  =4  G 
NNN  =  I N  C ” 1 

CG  nc  11=1*3 

XLC II ) =XC  G  00 (MM*  11  )) 

YLU  I  )  =  Y(NOO(MM*  II  )  ) 

ZLCI  l  )  =  ZCNC0CMM*I1  )) 

IOC  CCNTINUE 

>LC9)  =  “(  XL  C  1  )4  XL  C3XXLC5  )■•  XLC  /  ))/4.«(  XLC  2)4  XLC  4)4  XLC  6  ) 
1*XL( 3)  )/?. 

YLC9  )  =  -CYL (1  )*YLCi)*YLC5  )<YLC  Z))/4.*( YLC  ?)«YL(4)*YLC6) 
1« YLC6  )  )/2. 

Z  LC9  )  =  -C7LC  1  )♦  ZLC  3)4ZLC5  )4ZLC  Z  )  )/4.*(ZLC  2)<ZLC4  )42LC6) 
I  ♦  2L  C  2  )  )/  2. 

Il=rHCMM) 

IFCM SHELL. NE.O)  GO  TO  54 
CALL  SURVCCCTT) 

GC  TG  59 
54  CGNTINLE 

CC  36C  1=1 *N9 

VIC  I  )  =  ■). 

<2(1 )=0. 


36  C  V  2  (  I  )  -  T  T 
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59  CONTINUE 

NFLY=NHYER  IMH  ) 
rid  )  =- TI/2. 

00  J5C  1  =  i /  NPL  Y 
350  HCI )-H( I  )4THLC I#MM) 

OC  40C  LL=1»N^LY 
CC  5  K  =  I  /  o 
OQ  5  L = 1 r 6 
5  ZHX  »L  )=  ). 

■1T  =  mT  YPE(LL/HM  > 

IFClNC.CO.l  )  00  TO  6001 
3Ni2=bUH5TN(MM*LL) 

X 1=  A 1S(SN1 2) 

CALL  FONCT  (  000(1  )/ 000(51  >/ 000 ( 1 01 ) / ODOC  l5l)/rl*NN3/Xl 
1/  YP/  f PP/E5  ) 

3  KMT  )=  YP 
52<N  T) =YP 

6  CC 1  CONTINUE 

Al=WRANG(M T > 

ALl  =  Wr.'ANG(MT)*  3. 14159265  4/180. 

4  A  A ( l  )=  (H( LL  >*H (LLa 1)  )  /TT 
EL=£(HT> 

E1=E2(mT) 

VlL  =  ?f*  (MT  ) 

VTT  =  PS2(MT  ) 

3LT=G1(MT) 

5  I  =  5  2 ( “ T ) 

1F(CT.EQ.0. )ET=EL 

1F(GLT.E3.  j.  )3LT  =  EL /(?.*(  t.*\l  TL)> 

IFd/Tl.EO.Q.  )  /  T  T  =V  TL 
1F(3T.£0»0.)3T=ET/(2.*(1.<VTT)1 
Ud3=  5L  T 
313=51 

IF(  INC. ED. I )  GO  TO  302 
IFC  =  1 

CALL  YIELD (MM) 

302  CONTINUE 

CALL  OMAT(Q) 

OC  521  1  =  1  ,(■ 

CC  52 G  0  =  1/  6 

EEli  /  j)=EEU/J)4  3(I/J) 

52 C  CONTINUE 

DC  50  1=1,6 
CC  50  J  =  1 »  N  4  3 
i£(  I  / 0  >  =  ). 

50  3(I/J)=0. 

NNK  =  t 

IF(NELSP.NE.O)  NN«=2 
i/C  a  33  i(n  =  l/NNf( 

5  =  0. 

T  =0. 


in 
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IFCKK.EQ.l  )  GO  TO  601 
S  =  5 1  (NELSP  ) 

T=T 1CMELSP  ) 

SOI  CGNTIN  UiL 
7K  =  A4A<  l) 

CALL  SHAPE (S*r,7K*0CT J4C > 

OG  52  1=1*3 

l  f  CN  SHELL. ME  .0  >  GO  TO  53 
Ull-QiU  ) 

G21 =02  Cl  1 
J  21=0  3(1  ) 

0 12  =  “Vt  1  {  I  > 

G22=-»2CI  ) 

'J22  =  "‘«*i<I  ) 

GC  TO  57 
5  2  CONTINUE 
uii  =  o. 

U21 =1 « 

631= J* 

012=  1. 

G22= 0. 

U  22  =  0 • 

57  CCNMNUE 

call  eMATib*oo) 

52  CONTINUE 
CC  51  K  =  l  *  5 
OC  61  M=1*N45 
61  36<K*H)=Zrt*83CK*M> 

CC  5  L9=l*fi 
H  =  NJC(MN*L9  ) 

1<=5*CL9-1  ) 

CC  3  L  3  =  1*  5 
L  5  =  L  2  *L  3 
CCL5 ) =0. 

L4  =  I0(Ll*L  3) 

IFCL4.E0.0  )G0  TO  3 
CCL5 )  =  U(L4  ) 

3  CONTINUE 

if  (KiC.EQ.2)  GO  TO  630 
CC  600  K=l  *6 
01*1  =  0. 

SC  60 C  L  =  1  *  N  45 
u'wN  =  0UN*3(K*L)*0(L) 

0UM  =  0'JM*33  (K*L  )  *OCL  ) 

SCC  5IN(K)=0UM 

IFCNSrtLL.  Efi.O  ICAll  STRANSCSTL*STN»4L*0*0) 
IF(N5HELL.E0.0)G0  TO  610 
OC  62C  1=1*6 
62C  ilLC I >  =  jTN (  I  ) 

6  I C  CCNfINGE 

OG  S25  K  =1  *  6 
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QwM  =  0. 

CC  625  L  =  l  ,f: 

CLM  =  3LM«£.£(K.l)*$TL<L> 

625  SIGtfO-OUM 

CALL  STftANS(ST*SIG*AL»l»l) 

IFtlNCNAX. tQ.l >*RI T£( 6, ?000 H M>LL* ( 5  I  SC  I >*I  =  1 »6) 

GC  TO  6  33 
6  3  C  CCNTINUE 

OC  64C  K  -1  ,6 

0  lrt  =  0  • 

CC  64  C  L  =  1 >N45 
CUM  =  0Lf4*l(K»L)*0(L) 

0 IM  =  0UM'»  33  (*»L)*C(L) 

640  5  TN2 ( X  )  =  I)UM 

IrCNShELL.  ee.C  )  CALL  S  Tk  ANS(STL2*STN',»0.*G»C) 

I  ►  <N  SHELL*  E  T.  'j  )  30  TO  65  0 
CO  645  I=l»6 
645  3  TL2  (  I  )  =  STN2(I  ) 

65  0  CONTINUE 

00  oo5  K=l»6 
0lM=0. 

OC  655  L=i>6 

0L‘1=  3Lrf*EEtK,Ll*STL?fL  > 

555  SlS2(K)  =  QUri 

IFC  HCMAX.  EQ.l  )  WfU  TE  (6  #  2000  )  «M*LL*<  SIG2U  I  »I  =1*  5) 

63!  CONTINUE 

Ei«A-STL?n  ) 

allEXCMM^LL  )=5UHEX(NH»LL)*EXX 
S  >2  =  i TL< 1 > 

SLMSX2  MX,  LL)=SUiiSX2(MN#  LL)«5X2 

GSNl  2=2.*<  STLC2  )-STLC  l )  )*5N*C  S*STL(  4)*(CS  2-SN2  ) 

STM1  )  =  ST(  l  T 
5T*C! )=STC2) 

S1P(3)=ST(4) 

5  IGG  (  i  )  =  $I  G  (  l  ) 

S1GGC2  )  =  S  I  G(  2  ) 

3 IGG( 3 )  =  SI G  <  4  ) 

SUMS TN (MM,  LL)=$UMSTN< MM»LL)«0SN12 
CC  6  35  1=1  #3 

SLMS<  MM»LL  # 1  >  =  SUM5(.MM»LL,1>4STRCI> 

6  35  SUMS  1 G (MM*  LL  »1  )=  SUHS1GC  MM*LL*  I  >•»  SIGG(  1  > 

2  00C  FCRMAT  C2I5*18X*6E15.4  ) 
l  o3  2  FCPMAKIX,  IS.6C18.6) 

4CC  CONTINUE 
RETURN 
END 

SUBROUTINE  LOCAL(PRX*PRY,PFZ*  AA#NW*MM) 

CCMM0N/A/Q<9 ),XX<9  )*YY<?>,  ZZ(  S  )*V1<  ?>  »V2{  9)  */  3(  9) 

C  CMM0N/&/4HNU  53*5  >*PNN<  44,4  > 

CC*  MlN/M2/XX2(  9 ) »  Y  Y 2(  ■))  »7i2(9  )»Wl(?)*W2(9)»w3(9  ) 
ClfcNSlON  XL(9)»YL(9)#NFo0IN(4,3)*  XX1(9),XX3 
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US)#  YYIC9)#YY3(9  >#2Zl  (?  >,  2  2  3(9)#  PflX(l)#PRY 

2(1  )#PRZ(I  )  #  ZL(  9  )#  PLXC  IC)»°L  Y<  10)#»>LZ(  10)» 

3  KX( ?)#Wf( 9)#w?(9)#  ME LEM ( 3#  3)#NFP<9) 

4#AA(1  )#WW<  1  )  #PL1  (*?>#PL2(9>»PL3(  9) 

QATAN~LGM25#6»?#  4  #  0  #  3  #  3#2 

l#  1/ 

C  AT  4  NF-»OlN  /:>#7#1  #3#  6#«#2#4# 

i  7#1#3#52 

CAF4  NFP/l #2#S#6#9#3»7 #4  #5/ 

CC  U  T  =1 #  9 

PfcX( I >=0. 

PfiY(  I  ) - 0. 
f  FZ  (  I  )=}. 

11  CGN  T I N  oE 
CC  1C  1  =  1#  9 

X  >  3  (  I  )=XX2  (IHWl  (I  >*.5 
>>1(  I  ) =XX2(1 ) “ H 1 ( I >*.5 
n3iI)*YY2CI>«M2Cl  )*.5 
YU(I)  =  YY2(I)-*2(I)*.5 
2  2  3(1  ) =222(1 >«  *3(1  >*. S 
2  2  L  (  I  )  =2  22  (  I  )-*  3(1  )*••=< 

10  CONTINUE 

DG  12  I N  =1  #  4 
PN=PNN ( MM#  IN  ) 

IF(PN.G').  j  .  >G0  TO  It 
1  1=1 

12  =  4 

13  =  2 

DG  7  j=l#3 
K=NFP0IN(1 N# J) 

3  X  L  ( 11 )  =  XXt (K  ) 

Aid?  >  =X X?  (  K  ) 

Xl<  1  3 )  =  XX  3 ( *  ) 

YL  (1 1  )=  YY1  (K  ) 

YL(12)=YY2(K> 

Yl(K  3)  =  YY3(K) 

2UX1  >  =  ZZi  (K  ) 

2L(K2  WZZd) 

2L<K3)=27I(K) 

X I  =  1 1<1 
12  =12  *  1 
12  =  13*  l 

7  CONTINUE 
CG  1  *»  i  =  l#9 
IsNF’X 1 ) 

>)CI ) =  XL( K  ) 

YYCI )  =  YL(1  ) 

2  2(1 ) =?L(1  ) 

14  CONTINUE 

CALL  SURFAC(XX#  Y  Y#  ZZ#  V 1#  V  2#  V  3  #  l .  ) 

JC  121  1  =  1  #9 
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?LX( 1  )=0. 

PL  Y  (  I  )=,. 

PiZ< 1 ) =0. 

1 ? 1  CONTINUE 
AH4  =  C. 

CC  30 C  K«=!#2 
DC  3  0C  II =1*? 

:c  loc  j j= i  #  2 
S£-*MII) 

T 1  =  A  A  C JJ) 

Ztf^MCKX) 

*u  -*»<n ) 

HuJ  =  '4M(  JJ ) 

H«i(s  MhtKK) 

CALL  SHA°?(SS*rl»ZK»OETJAC» 

> >  L.L  1  =  wl  I*  WJJ*WKK*3£TJAC 
Ai»t  A  =  A  F<  t  A  f  XmUl  i 
DC  311  1=1#9 

'/  5  =  C  VI  (1  2 (I)**? 4  i  3(1  )**2)**.5 

w  i  =  v ;  c  I )  /  v  s 

i/*?=V2<l)/VS 
VY3=y3U  )/V s 
aLlU  )=PN* YV1*0C1) 

FL2(  I  )  =  PM*  YY2*'jCI  ) 

PL 3 1  1 )=PN* YV3*fl(  I) 

311  CONTINUE 

CC  312  1=1.9 

PL  X  (  I  )=PLX(  I  )♦  P  LI  C  I  )*  XM  UL 1 
PLY(I)=PLY(1 >4  PL  2(1 ) *  XMUL1 
PL(  (I)  =  PLZ  (I  >4PL3(I)*XHULl 

312  CCNTINUE 
300  CONTINUE 

pi.X  (  10=0. 

°LY( 10  >=0. 

PL7< 10  )=0. 

OC  321  J=  l  »  3 
«=NEPCIN< IN# J) 

Kl=N"LEM(i #J) 
f<2  =  N£L£M(  2#  J> 

K  3  =NC.L  £H(  3  #  J  ) 

1 F(K?.Efi.C  )  K 2  =  1 0 
3R  XI =PLX( K 1  ) 

PFX2=PLX(S?> 

FFX3=PLX(*  3  ) 

PfiX(<(l=PFX(*)<Pf<Xl«PPX2+PfcX3 
PFYl=PLY(Kl ) 

PFY2  =PLY(K? ) 

PPY3=PLY(K3) 

P?Y(  X  JspS  Y(!04PFri*PPY24PPY3 
P  F2 l  =  PLZ (K l ) 
aF/2  =PL  2  ( K  2  > 
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PB/3=PU(K3  ) 

°R7('(>  =  Pi’/  (X  )*PP  ZHPRZ  ?-*  ?9.7  5 
3°  1  CCNTlNUf 
ie  C  CN  T  I N  ut 
17  :CNT1NU£ 

RETURN 
E  M 

SU.ROlTlMt  3HAP2CS,T»7K  ,  0E  TJAC  ) 

CC9M0N/A/N  C  )),XL(9  )  ,YL<  7  >*2L(  9  Wl(9)  »  72  (  )  W3  <9  ) 
C1MEN310N  3I(a>*Tl(d),MS(9),MTc9),JA;i  3,  3) 

CAT  A  SI/ -l.,0.»i.,l.»l.*0.»-l 
DATA  T?/'al»,~i«,*‘la,'0*,l*,l«,ls,  0«  / 

?  I A ,N  »  N  S  »  4  T  ,  J  A  C 

CC  !  1  J  *1#»J 

DWl=l.«ii  (I  >*3 

01*  2  =  l.«Tl  (  I  >*T 

IFC1.£Q.2. 0F.1.ZG.6)  30  TO  13 

I  F ( I .£0.4. OF. I. EQ.fi)  GO  TO  17 

aw  $  =  S*  51 (  1  )4T*TICI  >-t. 

N  ( I ) =C  UM 1  *  DUH2*  0UM3/ A  * 

NS  C  I  )  =  ( OUM 1  *  GUM  2  *0UM2  *  0UM3)*Sl(T  >/m. 

N  T  C  I  >  =  <0W1*0U*1 2  ♦  DliHl* DOM  3) *T  1(1  )/4. 

GC  TO  11 
13  CUM4=l.-5**2 

.*4(1  J  *3W’* 00*4/2. 

N  S(  I  )  =  -  5*0  U*2 

Mr  I  >  =  T  1 C  I  )*OJM4/2. 

GC  TO  11 
17  CLM.5  =  i.-T**P 

NCI  )  =  CUMl*0UH5/2. 

N 5 (  I  )  =  0UM5  *  5  I{  T  )/2. 

N  1  (  1  )  =-T*L)UMl 
11  C  C  4  T i NU t 

N(9)=(l.-S**2)*(l.-T**2) 

NS(9)  =  -2.*S*(l.-T**2) 

M  (  7  )=~Z.*  T  *  (1  *  mS**2) 

OC  15  1=1,3 
3C  l 5  J  =  l»  3 
15  jfl  Zl  I  ,  j  >=0. 

CC  Zl  1=1,5 

JAC<  1,  1  >=JAC<1,  1  )*N5(  I  >*XL(  I  )«NS(H*7K*#1C1  >/2. 
JAC<l,2)=JAC(l,2)«NS(I)*YLCI)4NS{I)*ZK*V2(I)/2« 
J*C<  1  »l  )=J  *.C(1 ,3  MNSC  I  )»Zl(I  MN3(I  >*?«M  3(1  )  /  >. 

J  AC  (2,1  >=JAC(2,l)CMT(  1  )*XLCl)«NTCI  Wft*V  1(1  >/2. 

J AC( 2,2)= J 4CU, 2  >*NT( I >*YL(1I4NT<1 >*7K» V2(l >/2. 
JAC(2  ,  3  )  =  J  AC(2,  3  )4NT(  I  )*  2L(?MN  T(I  )  T(I  )  /  2. 

JAC(i,l)=JAC(3,l)4NCn*YtCI)/2. 

JAC<  3,2>=JACC3,2MN(I  >*V2(I)/2. 

JAC(  3,  3)  =  J4C(  1,  3)4NC!  >*<f  3<  I  >/2. 

n  cc  minus: 

J  £  T  J  AC  =J AC  C  1 , 1  )  *J.AC(2, 2)*JAC(  3,3) 


♦ JAC( 2,1) 
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1*JAC(3»2)*JAC(1»3)  4JAC(3*1>*JAC(1*2)*J4C(? 

2>2)  *JAC(1*3)*jAC(3*1)  *  JA  C  (2  *2  )  -JAC 

3(1*2 )*J*C<?*!  )*JAC(3»3)  -  JAG  <2*  3)*  JACC  3*  2) 

4*  JAC (1*1) 

*  AG 

ildSTUTINE  QMAT(Q) 

CC'MON/Y!  - L9/LL*rtT*EL*ET*r,LT*vTL*tfTr*G23*  G1  3*  ALl*  IPC 
CUSNSIUW  Q(S.S) 

OC  42,  1=1*6 
Ou  4  2C  J=i*f. 

420  C(I*J)=0. 

GS=CaSCALi  ) 

3f»  =  S INCAL1  3 
€32=03*05 
CS5=CS2*Co 
CS4  =  CS2*CS2 
3K2=$N*SN 
$M  =  SN2*SN 
5N4  =  SN  2  *  Sr*  2 

XMl./(EL*  (  1.-VTT**2)-2.*ET*VTL**2*  (1  .♦  rfTT  )) 

C11=EL**?* (l.-VTT**2)*XK 
C12=ET*EL*/TL*C l.*VTT)*XK 
Cl 3=C12 

C23  =  ET*(EL*VTT-*ET*VTL**2)*XK 
C22  =i T*CEL-ET*/TL**2)*XK 
C • 3=C22 

Cll=Cll-Cl 3* Cl 3/C33 

312=C12-C1  3*C 2  3 /C  3  3 

92?=C22-C2  3*C23/C33 

i fcb  =  3L  T 

t»  l  =9  12*2. *9  66 

U2  =  Qll-*322-A.*3Sb 

U3  =  Q1 1-Q12 -?.*3bb 

L4=0l2-Q224?.*Q6b 

'J5  =  Qli4a22-2.*Q  12-2.*Q66 

t(M  )  =  91 1 *CS4  4  2. *Ul *  SN2*C  S?4Q22*SN4 

3(1  *2)  -  li2*SN2*CS2-*31  ?*(5N4  4C  54  ) 

'3(',»2)=311*5N4  42.*U1*SN;,*CS24G22*CS4 
0(1*43  =  U  i*SN*CS34U4*SN3*CS 
w  (2  *4  1  =  J3*SN3*CS4U4*SN*CS3 
2(4*43=  Ui *  SN2  *  C  S24  06  o* ( SN4  4  C  54 3 
3(4*13=9(1*4) 

2(4*2  3-3(2*4) 

3(2*13=9(1  *  ?  3 

0(5*5  3=(G2  3*C5  24Gl3*SN2>*5./6. 

816*5  )  =  tr,23*SN2*Gl  3*CS2)  *5./6. 

C(5*6  )  =  (  (G2  3-.il  3  )*C  5*  SN>*5. /6. 

3(6*5) =0(5*6) 

*ETd«N 

END 


SUBROUTINE  SPLINECD»S»F#X#N#c:5) 

C IMF  NS  ION  3(S0),CC50),0<50),A(50),F<  3d)  »  SC  *0)  »XC  50) 

L  85  FORMA  I  (  ////»lX?bMQATA  FOP  CUBIC  SFLINE  FIT  / 

1  iX4  3P  NUM  3En  OF  SEGMENT  GI  GF  U 

Z  ) 

REAtH5*l3S>N*E*E5 
135  FCSM4IU5*  5X,2F10.0) 
hflTE(o»l?  7 >  N  •  £  *  E  '5 
187  FGFMAHI  5*  5X/2E2C.7) 

Nl=N4 t 

TE(6*2o  1  )  NT 

261  FOPMAKIOH  THESE  AftE*I2»7H  POINTS) 

‘Vf  ITEtfc*  3  30) 

330  FCrMAl  C57H  INPUT  MAG.  AN  0  INTERNAL  STARING  AT  l  ) 
hR J  TE  (  6»  340  ) 

24  C  FOPMAKbX*  4HFCI  )#17X#*rt><II  ) 

REAO  (5*  38  C  )(FCI  )#X(I)/T=i*Nt  ) 

2  8  C  FORM  A  H2FI0.0) 

MFJ TE( 5» 331 >  (FC  I )  »X<  I  )#  1  =1  »N1  ) 

3  e  t  FCPMAT < <2t 2C.7  )  ) 

OG  Zlf  I = 1 »  N 
S  (  T  )=X(H1  )“X(  I  ) 

217  CCNTINUE 

Ad )  =:. 

CCNl )=0. 

■j a  )  =  su  )/  3. 

u  oil  >  =  sc n  )  /I. 

00  54 C  1  =  1 » M 
C(J)*S(1)/A. 

UIU)--C(I) 
rC  =  1  *  l 

IF(l.EU.N)  GQ  TO  540 
3(N)=(S(IMS(K))/3* 

G(«)=(F(K4l )-F{K))/S(K)-(FCK)-FAK-l)>/SCI ) 

54  C  CCNTINUE 

C  <  I  )*(FC2)-FC1  n/S(l  )-E 
0(N1)=£5-(FCN1)-F(N>)/S(N) 

JG  IOC  I  =  ?  » N  1 
d(I)=-d(l)/A(l  ) 

C(I  )=-C(I > /A(J  ) 

0(1  >  =  -3(T)/AC  ) 

7  00  GENT  INGE 

CALL  SOL^E (0 »C *D»N1 *50 > 

RETURN 

END 

iLBRQUT  INF  FUNCT(0»  S,  F,  X,  Y*N  »Xi»yF»yPP»F5> 

C  I  ME  N  5  I  ON  0(1)#SCI)#FC1)»XC!) 

QG  2020  l = l »N 

lFCCl.liT.XCI*!))  GO  TO  2020 
U  =  0<  I  )*(x  (!♦!  >  -  X 1  )**?/(  6.*  SC  t)> 
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Y2  =  0{  Ml  >*  (Xl-XC  I  )  )**3/(6.*S(  1  )  ) 

ti  =  (F(I*l)/SCI)-CCl4i)*StIi/6.  »*(Xl-X<U) 

Y  W  r  (  I  )/  >  (  I  )-0  C  I )  *S(  1  >/6.  )*1  X<  1 4  1  )-x  l  ) 

Y=  Y1  4  I2-*  Y3  4Y4 

Y  1  -  “C  C (  I  )*  (  X  (  1+  l  )**X1  )**2)/<2.4SCl  )  ) 

Y2=C<  HI  )*  (  X 1  -  X  <  I  )  )**2/(  ?.*St  1  >> 

Y2=(F(  1*1  )-r(l  )  )/sfM 
YAs->(l)*tDC!^l)-l)(n»/f,  . 

Y  f  =  Y 1  ♦  Y2*  Y  J  4  Y4 

Y  l  =0  <  I >*<X< J 41 )-xl )/3CJ > 

Y«=au«u*cxi-xcn )/ >< ! ) 

YF»= Y1+Y2 

GC  TO  23S.) 

?C2  C  CONTINUE 

I F  C  X  L  .GT.  X(N*1))  YP=E5 
?  36C  CCNTIN  liti. 

RETURN 
£  A  3 

SUBROUTINE  SQLVE(fl*C»Q*Nl*ND) 

CIMENS10N  B<NO)*C(N01*D(NO> 

CC  l  2  4  j  1=1,  N 
N  =N  1  “  1 

s ( 1 4 1 j=aci >  *ac i ♦  i  >4C( n 

CCI41 >=9(I  )*C(I 41 J 
0<1  4 1  )=3Ci  )*DCl4l)4D(  I  ) 

I24i.  C  £N  T  I N  Ut 

C<Nl  )=v. 

OlNi  ) =DCNl  >/3CNl > 

cc  i  lao  j=  i»n 

I=.Ni-J 

XI  ) =  C 0< I >-C(I )*D< 1 41 )  )/5C I) 

13SC  continue 
return 

ENJ 

SUBROUTINE  Y IEL  0  (MM  ) 

COMMON /SET2 /INC* NFLAG*YL0XC4)* TIDY <  4>  *YLD$<  4)  *YL3XXC4> 

1  *  YL0YY(4)»SoNS<100»8»3)*K0jiNTC100)* 

2  XOUNTFUOO  )*NFLA1LU  00*  3  )  ,  N  Y  I  ELXl  00  *  3  > 

3*SUM$IGC  IX  #3*  3)  *  SUMS  TN  (l  CO*  2  0*  SUM£X<  100*30) 

4  *SUNSX2( 1 00*20 ) 

CCMMQN/YIELO/LL*MT*EL* ET*GlT*vTL*Y  TT*  G2  3*  G1  3*  ALl»IPC 
CCNMON /P/NS  HELL *N9*NFDG EL*  XOUNTY*  XOUN 
0  A  T  A  0,Nt*T«Q*3LAN/6HYEILD**6HFAIL***lH  / 

NNN=INC-I 
SS=YLCS(HT  ) 

3  I  T  =  Y  L  i)  S<  rt  T  ) 

Fl-1./ YL0X(MT)-1,/Y10XX<  NT) 

F2=l./ YL0YCMT)-1*/YL0YY(MT) 

F2=lwSS-l./STT 

rli  =  i./(YLOX(Mf  )*YLOXX(MT>) 

F22  =  l./CYL0YCNT)4VL0YYC»»T)) 


- pH 
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f 35  =  1 ./(SS*ST7  ) 

F 1 2=  9 • 

$l  =  S’Jl*StMM  /IL*  1  ) 

S2  =  5UMS(MM#LL*2  ) 

S12=5DMS(MM»LL#  3  ) 
i  liiX  =  SUMSl  G(MM»LL* 1  ) 

3  1GY  =SUMSI  G(MM,  LL#  2) 

SIGXY  =  $:J1>1GC>!M#LL*  3) 

E>X  =  jlV'.Ki  >M»Ll  ) 
iX2  =  SUM6X2(MM>l.L  ) 

IF<NFHS.ttt.2>  GO  TO  3*5 

°C=(  51/YLO  X  ( *1 T  )  ) **  2« (  S  2/YLO Y(  MT  )  ) **  2  -  5 ! *S 2/< YL  OX 

KMT  )  *  YJ.OXC  MT  >>•»  C  S  12/S  S)  **2 
uC  TO  312 

3C9  PC=F1*$HF?*S2*F  3*512  * F 1 1  *5  1* * 2«F 2? *S 2**2 

1  *  F  3 3*Si2**2*2.*F12*Sl*  S? 

312  CCNTINUE 
T  A '  i  =  3  c  A  M 

IFCPC.LT.l )  GO  TO  301 
I F  <  S  l  •  GT  •  >  •  )  GO  TO  5  j? 
x>X=- YL9ax(*T> 

IF( Sl.3T.XXX)  GO  TO  305 
3C  TO  300 

3C  /  IFC  51  .LT.YLDXCMT  ))  GO  TO  305 
3ca  EL  =  nc. 
c  T  =  1 9  0  • 

GIT=1C0« 

YTL=). 

YTT  =  ). 

12  3=1  CO. 

G1 3=1  CO. 

T  A  1=  T  MO 
GC  TO  3  0s, 

3C5  * T  =  ICO. 

GLT=lGO. 

Y  1L*0. 

Y  T  T  =  3  . 

553=100. 

T4G=0NE 

306  flOiJMlY  =  l 
301  CCM?  1  NU£ 

IFOPC.E3. 1  )WHIT£{6»5?9)  NNN»MM*TAG*LL*3C 

1  #  S  T IX #  5  l  GY  »SIGXY*Sl»52#S12#SX2*FXX 
55  0  FC9Y4TC1X»?I3*  Ab»l2*5XrP5.2#3tU.M 
»  £TU;«.\ 

END 

5 18  POUT  INE  INLQ  AO  (  X,  P  Y  *P  Z  *P  Ph  p  PU  *  N£  0  Ge  L  ) 

CCMM0N/inA9GH/L0*NEQ»NL0  AO #N£  L* NEIP L>  NNP 
DIMENSION  h  C302KPXCA4  ),p>C4«,  ),PZ(  4  4)  #PM{  44»4)  #PNC44) 
1  /  PU<44 * 3  )  *NMOE(  4) 

CCMMON/Hl  /  X  ( 1 G  3  )#  Y  Cl  63  )*  7(1  6  3  ) *  I 9( 1  63*5) 
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F~~ 


CAT  A  NS  1 0-:/ 3, A  .1,2/ 
mH  1  TE(a,tOf>!)) 

SC  T  3  (l,2*3,A,5*9,9,9»y)*LD 
A  CCNT1NUE 

mf> i  Tetb.n  si  > 

Oil  10  1=1, NE3 

l'J  S(1)=C. 

OC  2  2C  L=1  ,  NLG  A  0 
P  c  A  Q  (  5 , 1 00  0  >  ND  »  101  PN,rLnAO 
Ac  IT  £(6,10  Z  C>NO.  ID  I  Sn,  FLO  A3 
J I =1  ")(N0, I  0  I  PN  ) 

Iff  I  I  >220,  '>?0,  ?  AO 

2a;  fui  >=scii  mfload 

220  CONTINUE 
b  i  T  *J  rt  M 

1  CONTINUE 

OC  3?C  M  =1  , NEL 
Wn  IT  2(6,10/0 
*X(N  )=0. 

P  Y  ( ,-n  =  0 , 

PZC*O  =  0. 

39*  CENT  IMJE 
*N=0 

00  A  2C  1  =  1, NELPL 

IFCKN.E'J.j  )F5 *0(5.1 001  )L.PXCL).PY(L>*PZCL)  »KN 
IFCitN.EQ.O)  GO  TO  A10 
PHI)  =;»X(  t  ) 

FYCI  >  =  ?Y(i  ) 

Pit  I >=PZ(1  ) 

4lC  t»MTE(6,l0  01)L»?X(L>,PY(L)»P/(L).KN 
42 C  CONTINUE 
PETUPN 
5  CONTINUE 

W  F  I  T  E  (  6, 10*0) 

CC  3  3  0  N=1  ,NEL 
OQ  3  3 C  L  =1 , A 
32C  PMtM,L  )=0. 

OC  3ZC  H=l, NEOGEL 
PEA  Ot  5,1)1  0  )L.  TSIOE.PNLD 
uni  TEC 6, 101 5  )L, I  SI OE,PMLO 
NS  1 3=N j IDE  ( ISIDF  ) 

PM  t  L»  N  SI  D  )  =PMLD 
3/C  CONTINUE 
SETUPS 

2  CONTINUE 

WFI TE(6»10?P) 

*N  =  3 

OC  500  M  =1  ,NEL 

IFCKN.EO.O  )PFA0C5,10  30)L,PNCL >,KN 
IF  (KN.EQ. 0  }  GO  fC  450 
PNMMPNC1  ) 


r 
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450  hRlTG(6*l: 30)  M.PNLM)* KN 
3CC  CCN  T I N  LE 
P  £  T  ■  J  \ 

J  continue 

4F  1T.G(6»>000) 

CC  5  l  C  1=1  *NELPL 
«E*3('»rKm  L*<PtKL»J)»J  =  l*P) 

WHT"(«*10«i0>  L» (PUCL*J)*J=l*c> 

310  CONTINUE 
9  CG.NUNUE 
ScTURN 

t  C 1  5  FGt’rtA1(lX*I5,If*5X>F12.4) 

1CCC  fGPMATMO*  5X,3F10.4) 
l  30  l  F  Cf*M  A  T  ( 15*  5  X  *  3F1  C. 4 *  1  5  ) 

1010  FC~MAT(2I5*F10.4  ) 

1  C 3 C  FGftM4T(I5#5X/Fl0.4*I5) 

1C4C  FCi<MAT<l5, /,  8F10.4) 

1061  FCftMAlC  IX*  20HCONC  €N  TP  EO  LQAO  //// 

1  IX  ZiH  NODE  0 1  f  £  C  T 1  ON  LOAD  7 

.?  IX  32hNlH3Gfi  MAGNITUDE  //  ) 

1070  FCPMATt  IX  ZuHOISTPIOUTEO  LQ  40  7777 

1  IX  41HELEMENT  PX  PY  PZ  XN  7 

2  IX  10H  NO.  > 

FCPHATt  IX  20H  G  DOE  LOAD  77/7  IX 

1  40KELE.AGNT  SIDE  LQAQ  7 

2  IX  ACM  NO.  M AGNIT UOE  7 

3  ) 

l-)«0  FCrtMAH  IX  11HN0PHAL  LoAO  77 

1//  IX30HO.EMGNT  PN  XN 

2  /  IX  SH  NO.  7  } 

2000  FC.NMATC  IX  35HN0NU.N  I  FQfcM  QISTM3UTED  LOAO 

1  /  lx  «*8HELFm£NT  PRESSING  IMTFNSITY  AT 

2  LOCAL  COOK.  7  IX  45H  NG.  1234 

3  5  b  7  8  7  ) 

1C5C  FCPHAKIHl  IX  27HL  OAO  TYPE  DATA 

1  77/7  ) 

E*J 

5L1R0  L  TI N£  T NP* T C X  * Y * 7 * NNP I 
C  IMG  NS  I  ON  1P*CC1  >»XCif>M»VClSlW<165) 

C  AT  A  lpKC/l HF / 

RA0=ATANC1. 0>/(45.0) 

NOLO  =C 

rffl  TF (&,?.»  30  ) 

1C  continue 

READ  (5*1300)1 T*N*X(N)*Y{N)*Z<N>*XN 
1000  FCuMATC  Al,  14, 5X*  3F1  0.  0,1  5) 

nM  TE(6*?0  02)IT*N*X(N  )* Y(N >*Z (N)*KN 
2CCZ  FCPMAT<IX*A1*I4*3F12.4*15  ) 

IF(NHL  O.FO. 0  )  GO  TO  00 
I  F(iLN.E9.0  )  GO  TO  50 
NU  =  (N-NOLO  >/KN 
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r~ 


NtHN=NLM“t 

! F (NUMN.LT . *  )  GO  TO  50 
XMJM  =NUM 

CX=  T  X  (  N  )-X  (NOLO)  )/  XNUM 
OY=(Y(N)-Y (NOLO )  ) / XNUM 
0Z=(Z(N)-7 (NOLO )  )/XNUM 
A  =  NfJL  i. 

DC  TO  J=1,NUMN 

ah^-k 

A  =A  ♦  K  N 

X(K)=X(K1)*DX 
Y(K  >=Y (KKMOY 
200  =2<XM)-»0Z 
3C  CCNT1.NUE 
50  NCLO  =N 

IFUT.NE.IPRCC1  )  >G0  TO  SO 

IF(KN.£<).0  )G0  TO  70 

CC  29  J  =  1 *  NUMN 

k=n-J*kn 

CLM  =  ZOO*F  AO 

Z(K)  =  Y(<O*C0S(QUM  ) 

Y(K)=Y(IO*S!N(OUM) 

2S  CENT INUE 
7C  CENTINUE 

OLM=Z(NOLO  l*R40 
Z  (NOLC  )=Y(NOLO  )*COS(OUH) 

Y  (NOLC  )  =  Y(  NOLO  >*  SIN  (OUM  > 

6 C  NGLD=N 

IF(N.NE.NNP)GO  TO  1C 
RETURN 

?0?C  FORMATE /// /, 1 X  22HNOOAL  POINT  INPUT  OATA  /  IX 

1,43H  NOOE  NQJAL  POINT  COORDINATES  / 

2  IX  4 3  HN'JM  8FR  X  Y  Z  AN/ 

3  > 

END 

3 COROUTINE  INIOC  IO,NNP,NEQ) 

DIMENSION  10(163,5) 

NOLD-C 

NFITE(6,?0  35  ) 

15  CONTINUE 

*£40(5,1005)  N,(ID(N,I ),I=i,5 >,*N 
1005  FCRM4K7I5) 

l»FIT£(6,?)C6  )N»  (10  (N,l),  1=1,5  ),«N 
20C€  FORMA  1(1 X, 715) 

IF (NOLO. €3. P)GO  TO  55 
3C  20  1=1,5 

IF(ID(N,I).EQ.O.ANO.IO(NOLO,I  )«LT>0)  ID(N 

1, 1)= IC(NOLO.I) 

20  CCNTIAUE 

IFCXN.E3.0)  GO  TO  55 
NUM=(N-NOLO)/KN 
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NLMN=NUM“1 

I  Ft  NU**  N«  LT  •  1  )  GQ  TO  55 
K=NQL0 

CC  35  J=l*NU*rt 

K»<=* 

K  =  K«  XN 
CC  55  1  =  1*5 
ICC*.  I  >  =  1 .3  (KK*  I  ) 

If  (I’3(K*I  ).GT.  I  )10<K*  I  )  =  lD(Kt#i  >4RM 

35  CONTINUE 
55  NCL3=K 

IFCN.NE.NNP)  53  TO  15 
Wnl  T  E  (  6*  20  4  G  ) 

WF 1 T E( 6*201 0 ><N*  C 1 OtN * l > *  I  =  1  *  5 > * N=1 *NNP) 

2  )lvi  FCMAT  (515  ) 

N  E  i  =  0 

2C  52  N=1*NNP 
CC  52  1=1*5 
IC(N*I  )=I45S(I0(**1  )) 
lFtiO(N*i)*l) 5 7*53*59 

5 1  AEC  =  'JE34l 
ID(N#I  )=NZO 
GC  TO  62 

53  IC(N*I  )=) 

GC  TO  5? 

59  U(M  »  1  )=-lD(N*  1  ) 

52  CCMINUt 

20  35  FGMAK///.1X  20H1NPUT  10  C  JOE  S 

1  /IX  40H  NGOc  90UN0ART  CONDITION  COOES 

2  /IX  4UH-N  LHBEP  X  Y  l 

3  ALPHA  9 A  T  A  XN  /  } 

2C4C  FC&MA  i(//i  X.40HGENEMTFO  10  COOES 

1  /  IX  4  OH  NO  Of  9  0  UN  0  49  Y  CGNOIUON  COOES 

2  /IX  4 ) HNUM3ER  X  T  Z  ALPHA  BATA 

3  /) 
fi£T JSN 

EN  J 

SUBROUTINE  SUFt/ECtTHICX) 

CCMM  3fi/H2/XL(  ?)*  YL<9>  »ZLC9)*y  1(  9)*V  2<  9)  *  y  3t  9)  *JAC<  3*3) 
l*  N(9  ),NS(9)*NT(9)*Ui(9>*U2C9>*U3(9>*Wl  (9),  M2  (9) 

2**3C?  ) 

DIMENSION  SI  (  4  )  *  TI  (  8) 

REAL  N.NS.NT.JAC 

?EAL  L 1  *L2  *  L  3  *  M  l  *M2*M  3*  N1 *N2*  N  3 

DATA  SI /-l  .*0.# !.*!.* 1**0. *-l .#-!•/ 

DATA  U/-U  *-l.*-l.*0.*l.  *1.  *1.*0./ 

CC  30  11  =  1*? 

$*51(11) 

T*T ! ( 1 1  ) 

CC  ll  1=1*3 

cuni=i  *4sm)*s 


0UH2=i.*Tl  <  I  I*  T 

1FC1-EQ-2-  Of  .  !.  EQ-  6)  GO  TO  13 

IFC l-EQ-4. 0F.1.E9.8)  GO  TO  1 7 

dlM3  =  S*  SI  (  I  )♦  T*TI(I)-1. 

MI  )  =  CUM1>DUH2*0UM  3/4- 

N  S  (  1  )~  {  OUii  1  *  OUM  240UH?*9UM3)*$T(I  >/4 

NIC  I  )-(D,JMi*OUM2*OUMl*OUM3)*TICT  )/4 

•SC  TJ  li 

13  CUM4  =  1 . ~  5**2 

NCI  )  =0Ui“l2*l.)lM4/2. 

N  SC  I  >  =  -  5*0 GH2 

NIC  I  )  *  T  I C  I  >*i)U*4/2. 

GC  TO  11 

17  04M5=1--  T  *  *  2 

NCI  >=CUMl*0UM5/2. 

NSC  I ) * 0UM5  *  5 T C I  )/2. 

N  1 C  I  )  =“  T*0UM1 

11  continue 
XXI*')- 
ym  =o- 
1X1 *0. 

X  E  T  =  0  • 

YcT=9. 

7  E  T  *  0  • 

GC  20  1  =  1,8 

fl>l  =  XXl4NS(l)*XLCI> 

r>i=y>i4NSCi ) * yl < i  > 

ZXI=ZXI4N6  CT)*ZLCI  ) 

XET  =  XET4NTC  M*XLCI  ) 

YET  =  YET4NTCT)*YLCI  ) 

Z£T*ZE  T  +  NT  C  I  )*ZLC!  ) 

0  CENT INUE 

QL  =  t  XXI*XXI*YXI*YX1#ZXI*ZXI  )**-5 
L1=XX!/DL 
M 1 =YX 1/ DL  . 

N 1  =  7  XI /DL 
J 1  <  I  I  )  =  L1 
02 C  I  I )=H1 
U3CI  X  )=N1 

JL  =  (XET*XcT*YET*YET*ZET*ZET  )**-5 

L2=XET/0L 

M  2  =  Y 1 1  /  OL 

N2  =  2ET  /DL 

W1C1 I  )  =  L  2 

«i2C  I 1  )=f12 

k  2 ( X ! ) =N2 

L  3=M1 *  N2“M  2*Nl 

M3=L?*NI-N?*L1 

N!=Ll*M2-«l*L2 

V  i C 1  I  >  =  L  3*  TH ICK 

421  II)*rt3*THICK 
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V3  ( 1 1  >  =  N3*THICK 
3C  CGNTINUE 
RETURN 
END 

SUBROUTINE  ETRANCE.0E.0) 

CC*S*3N/H2/XL<*>.  YL(*)  »ZL(9).Y  1(9>»V2(  9>*V  3(  9)  , 
l  JACC 3.3>.N(9).NS(9).NT{9) 

CCMN  0N/0C/L1  .L2#L3#M1.M?.M3.nI»N2*N3 

0HENS1G*  E(5.b  ).D£(6.6  ),  0(6.6)#  3  TEC  6.  6 )  .OE TE(  6  »M 
1.IE(6*6).  E  T£ { S .  6  ) 

SEAL  N.NS.NT.JAC 

REAL  LI.L2  »U.*1.'«2>"  3.N1.N?»N3 

)  >  T  =  j  . 

Y  X I  =  3  • 

/XI  -  ) . 

>  e  r  —  3 . 

Y£T  =  0. 
z  t  r  =  i . 

ZL  2  3  1  =  1.  *» 

XXl  =X*I«NS(I  )*XL(I  ) 

YXI  =  YX  UNS  (  1)*YL(I  ) 

Z>I=Z>  HNS  (I  >*7L(I  ) 

Xtl =  XE  T*NT (I )*XL(I  ) 

YtT  =  YIT4N'T  (  I  )*  YL  (  I  ) 

2ET  =  7ET*Nr  (n*7LCI) 

2C  CONTINUE 

CL  =  (XX1*XX:  +  YX1#YXI4/XJ  *ZXI)**.5 

L’=YX1/0L 

Ml=rxi/DL 

N1=/XI/0L 

CL  =  (XET*XET«  YET*YET«7ET*ZGT)**.5 

LE=XE 1/DL 

*2  =  YET/[)L 

N2  =  Z  ET/OL 

L3=M1*N2-«2*NI 

M2=L2*NI-N2*L1 

N  3  =  L l *  M  ?-M  t  *L2 

L2  =**  3  *  N!  “N  5*11 

(2=L1*N3-L  >Y1 

N2=L  5  *  ill  “  3  3*L1 

CALL  TMAX(  TE.l  .  #  2.  ) 

JC  IOC  1=1.6 
3C  IOC  J= 1 . 6 
El£« 1»J>=). 

Dc TE (  I .  J )=  i. 

C f £ (  I.J)=0. 

OC  IOC  M  =1  .6 

£1EC  I»J>*£TECI»J  >«EC!»*0*TE<Jl»'l) 

OETtC  UJ)=0CT?C1»  JMOEC I.HWflM.  J> 

CTE(I.J)=  3  *E  (  I  .  J  !♦  OU»*n*Tt(H.j) 

ICC  CGUTINUf 


W  CO  23C  1  =  1*6 

*  CC  2-VJ  j=t*ft 

"  £U*J>=0. 

QEt 1  * J  )=0. 

CCI*J>=). 

CO  2  33  M=1  *6 

Ed  *  J>  =  EU  •J>*TE<M*I>*ETECM*J  > 

C£<  I*J  >=0£  U*J  )«T£(M*  I  MDETEf  **J> 

C{I*J)  =  CCl*J>«TECM*I>*OTEtM*J) 

2CC  CONTINUE 
RETURN 
£NO 

SUBROUTINE  STR  ANS<  TT» S1G*  AL  *N  N *  KKK) 
CCMMGN/H2/ X  (9  >*  Y  (9  ),Z  (S>  ),  VI  C9  >*  V2C  9  ),  j(  9)  * 
1  „AC(3*i)*N(9)*NS<9)*NT<9) 

1776  CONTINUE 

COMMON  /OC/L 1  *L,*L3*Ml*M2*M3*Nl*N2*N3 
DIMENSION  T<6*S)*TT(6)*SISt6> 

REAL  N #NS#  N T*  j AC 

PEmL  L1*M1*N1*  L2*M?*N2*  L3*M3*NS 

If tXK.EG.OJGO  TO  2 
P  1  =  3.  H159265A 
a*AL*PI/l*I9. 

L1=C0S(8 ) 

K1=STN(3> 

ni  =  :. 

LC=-SIN(3) 

M2=COS(3) 

N2=0  . 

13  =  3. 

M  3  =  0  . 

N3  =  l. 

G  C  T  J  3 
2  CONTINUE 
Q  >  =  3  • 

0  =  0. 

Zi=  o. 

CC  13  1=1*8 
3X=  OX^NSU  )*X(  I  ) 

C»  =  Oy-»NSU  )*Y(I) 

13  Q  2  =  0  Z«  N  $(  1  )*?(!) 

OL  =  <  D  X  *DX*  QY*0Y*QZ*DZ)**.5 

oi-nx/oL 
M 1  =  0  Y / OL 
N1  =  D  Z/OL 

cx  =  o. 

0=0. 

02  =  0. 

CC  20  1*1*6 
C>  =  0X*NT(1  )*X(I  ) 

0  Y  =0  Y  ♦  M  T  C I  )  * Y(  I  ) 


Oi  =  T?4*T(  1  )*?(!  ) 

CL=(  3X*DXOy*.>  Y<OZ*PZ)  **.'5 

Ic  =  OX/QL 

Ac  =or /OL 

.\!  2^=  07  ✓  L>L 

L3=M1  *N2-M?*M1 

M  3  =L  •  Nl "N  ? *L  1 

N3  =  Ll*M2-d*L? 

12 =  M3 *N1 -N  3*MI 
M2  =  LI*Nj“L  3*N1 
N2  =  L  3*Ml”M  3*cl 
3  CONTINUE 

I FCKKK.EQ. 0 >  <30  TO  5 
:  =  .?.  ' 
c=t . 

GO  TO  & 

CGNT INUE 
C  =  1  . 

0=2. 

CONTINUE 

CALL  T  M  A  X  (  T  /  C  /  0  ) 

GO  1  I=l/o 

rid  )=o. 

OC  l  J  =  l/*» 

TUT  )  =  TT( I  >«T(!/J>*$lu(J> 

1  CONTINUE 
RETURN 
END 

SlididUTINE  TMAX(T/C/0> 

DIMENSION  T ( 6»  6  ) 

CGMM0N/DC/L1/L2/L3/M1/M2/M3 /N1/N2/N3 

REAL  LI.L2 /L3/M1 /M2/M3/N1/N2/ *3 

TCI  /  1  )  =Lt**2 

T  (  1  »  2 ) =Ml *  *  2 

TCI  /3)=Nl**2 

T(l#4)=C*Ll*Ml 

T(lz5)*CMl*Nl 

TCI »6 )=C*N1 *L1 

TC2/1) =L  ?*  *  ? 

TC2/2)=M2**2 

TC2/3)=N?**? 

T(2/4  l*C*L2*M2 
T(2»5)=C*M2*N2 
T  C2 z6 )  =  C*N 2* L2 
1(3/1 )=L3**? 

1(3/  ?)  =  M3**2 
T(3/3j  =  N3**2 
1(3/4  >  =  C*L3*M3 
I  ( 3/  5 ) =C*M 3*N  3 
r ( 3/6)*C*N  3*L 3 
TC4  .1  )=L1  *  L,*’J 
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T  (4  .?  )=tfl*M?*0 
K4,  3)=Nt*N?*l) 

F<4#4  )=Ll*4?*L?*Ml 
T(4#S  )=M1*N?«M2*N1 
T (4  »6  )  =N1 *L?*N2*L1 
T (5. 1  )  =  L2*L  3*0 
T<5  ,Z  )=M?*M  3*0 
T(5,(  >  =  X’*N3*0 
1(5 »  4  > =12**  3*L  5*M2 
/  T ( 5* 5 )=H2*N 3*M3 *N2 
r(b#6  )=N2*LMN3*L2 
i(6*  1  )  =1.3*  LI *Q 
?>=M  3*  M  1  *0 
F(b#  3  )=N3*N1*U 
T(o.4)=L3**tl-»Ll*H3 
F  (a.5)  =  M3*Nl-*Ml*N3 
F(6»6  >  =  N3*L1*NI*L3 
f<  £  T  l*  S  \ 

ENG 

s ueo  11* t  i N£  m  a  r  (3*88) 

GI-rINSIQN  d(6.45),33{6*45)#KKlU5) 

CGHM0N/K/K1  .K2,K3»K4#K5 

CCM*nN/H?2/TT,ull#U2l,U31,  UI?*U22*U3?  *1 
CCMir]N/H?/XL(9),  YL  <  9  )  ,  7  L  <  9  )  *  / 1  <  9  >  ,  V  2(  9>,V  3{  «)  * 

1  J A C ( 3  »  5)»NC9|#N3(9)*NT(9) 

r’EAL  jac.n^ns.nt 
OC  51  J=l*  5 
51  K KK ( J ) -5*( I  - 1  )*  J 
lU  =  KK*(l  ) 

r.2«K«KC2) 

<3  *  =  KK^C  3) 

*4=KK!((4  } 

«5-K«!U5) 

Al-JAC<l»l)*NS(l)*JACCl»2)*NTi!) 
31=J*C(2»l)*NS(IHJAC(2,  2)*NF(I) 

C 1=JAC<  3*1  )*NSC! H  JAC( 3»  2)*NT (1 ) 

3  ( 1  *  K 1  )  =  4  i 

3(l/K4)=JACCl»3)*M(l)*Ull*TT/2. 

U1  '\i  )=JACC1»  3  )*N(l)*U12*TT/2. 

3  C 2 #  K  2  >=3l 

3  (2#  <4  )=JAC  (2*  3)*NCi  >  *U21*TF/2. 

3(2  j  )=JAC(2»3  )*N(I)*'J2  2*TT/2. 

3C3/K2  ) -C 1 

3  (  3»i<4  >*  JAC  (  3*3)*NC1  >  *U31*TT/2. 

2  0  ^5  >=JAC<3#3  >*NU)*U32*TT/2. 

3(4*IU  )=0l 
6 (4#  K  2  )=  A  l 

3 (4 . K4  )  =  JA  C (2*  3 >*N( 1 )*  Ut  l*TT/2.*  JAC  <1  *  3) *N( I ) *U2t «TT 

1/2. 

3(4*K5  >  =  JAC  <2*  3  )*N(I>*(Ji2*TT/2.*JAC<l  *3  >*N{  I  >*(,22*  TT 

1/2. 
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3(5,K2)=Cl 
3 (5*  )  =  R1 

3<5,K4)=J\C(  S, 3 ) *N(  I)*'<*21*TT/2.*JAC(2*3)*N(  !)*U3t*TT 

1/2. 

d(i,K5>=JAC<3,3  >*N(I >*U2  2*  TT/2.* JAC  <?,S>  *N( l ) *U32*TT 

1/?. 

<3  t  c  » f  1)  =C  l 
3  (fc  *X  l  )-  A* 

3  (Hi.  *4  )  =J4C  (  3,3)  *N(  !  )*0l  1  *TT/2  .♦  JAC(!  ,  3)  *NC  I  >*1,  31  *TT 

1/2. 

3(b»*<5  )=JAC(S*S  )*M( 1>*U1 2* TT/2 .4  JAC<1  »3>*N(l)*U32*TT 

1/2. 

3SU*K4  )  -  A  i.  *  J  i  1  *  TT/2. 

3fc  (1  ,K5  )=  4i*Ul  Z*  TT/2. 

3<3(’,K4)=31*U?l*TT/2. 

36<2»K5>=3  1*U2  2*TT/2. 

3  £  (  3,A4  )=C  l*'J31*TT/2. 

3  £  (  3  ,A5)=Cl*U32*TT/2. 

3B(4,X4 )=3 l*UU*TT/2.4  A1*U21*  TT/2. 

33<**,ri5)  =  ?l  *U  l  2  * TT /  ? « •*  A  1  *L»22 *  T  T / 2  • 

3E(5,X4 )  =  Cl*U2t * TT/2.  431  *U3l* TT/2. 

3c(5»K5)=Ci*U2  2*TT/2.43l*U32*TT/2. 

33(5*K4)=C  l*oil*TT  /2.4Al*Li31*TT/Z» 
3£(6,ft5)=Cl*U12*TT/2.4 4l*U32*TT/2. 

PET'JRN 

-NO 

31 3300  TIN1:  *FSH3(XX,YY,  ZZ,NQD  *>Nt3*NFL,.ND,  TITLE  > 
GIHEN3ICN  TJTLE:U3>,XP(ICC  ),YP(10„),XP3( 1 3) »YRG< 1 3 >,N 
l  ( 1 2 ) , NQN( 1 2) 

01PEN3I3N  N  N  (  2 1  ,21  ),  YC  (  2  l,  2  i )  *  XCC  21  ,2  l  >  *  NNP  3<  20*4,21) 

l*«iT(.?j,4  ) 

J  IMF  N  SION  L5(  d)  ,  NF  (  40  0)  ,XF(40C),  Yd  (  400),  NP ( * > , IC3MP ( 4 
1,4) 

C  IMBNSIQN  XX (i  ), YY(1 ),N30CN0, £ ) 

3  HENS  JON  2P(100),ZPG< 1 3>,Zf( 40fi),Z2C 1>,ZC<21 ,21) 

REAL  N 

3  A  T  A  I  COMP /- 1,1 , 1,-1, 1,-1, -l,  1,1, -1  ,-l, 1, -1,1 ,1,-1/ 

C  AT  4  N3W/(V,N3/0/ 

CAM  Z,U,Y.W/0.,C.,0.,0./ 

RErflNC  1 
NE«3  =  0 
NtL  =  0 
NEN  =  « 

»MTEC5,1M  (TI  TLEC  I)  ,  1  =  1,11  ) 

17  FCHATCiX,  MESH  GENERATION  FJfl*llA6,/> 

READ  (5,1)  iNrcG,  INBP,  N3N 
l  FGPMATC4I5  ) 

PE  A9 ( 5 , 3)( X  p (1 )#I=l*IN9D) 

«EAO(5,3)(  YP(I  ),  I  =1  ,  IN  3P  ) 

P£AD(  5, 3  >( ZP ( I  ),!  =  !,! NOP  ) 

I F  <  V3  N  • E  3.  3  )  N0N=  ft 
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3  FQfiMAT<aF10.5) 

CC  2  J=l# I NPG 

2  REA’)C5#3)  N9G#<  JT(NPG»J  )»J  =  1  #4) 

6  FORM  AT ( 51 5  ) 

XH  TE(b#  3b  ) 

36  FCF1ATUHO ////l  X#"GL08AL  COUP  01 NATE  S  "// lX  »  "NU  H8ER  X 

i  co 3 pc  y  coopo  "> 

Xfirr:(b#5:  )  (1  #  XP(I  )#  YP<  1  )#  l>  (  1  >#J  =1#  IMBP  ) 

3 C  FCPH \ 1  (?X#I 3»7X#F7.2#5X#F7.?#5X»1  7 .2) 

XPITE (6#  ’1  ) 

21  F  C  nN A  T  (//1X#17HC0NNECT1YI  TV  0  A TA/1 X #4 lH REGI ON  SlD-I 

11  2  3  4  ) 

CG  2*>  I  =  1#1N"G 

2  6  W  FI  TF  («»»??)  I »  C  J  T  (  I  #  J  )  #  J  =  1  #  4  ) 

22  FCF-IA  1(2X»  I  3#14X»4(  I2#5x  )> 

0G16KK=l#I 

FcA0(5#4)  NPG#NROWS#NCQL#CNON(I >#I=1#N9N) 

4  FCFMAT(1515 ) 

up  I  TEC  6#  )  MPG»\K0v*6  »NCOL#<NCN<  I  )#  1=  l# NON) 

Id  FCPMA  T  C 1  HO  ✓  ✓  /  1 X  #  1 2  H  *  *  *  REGION  #I*#6H  **** / /1QX# I ?#5rt 

1  FQMS»l  )X»  12*4H  C3LU"4NS//lO  X#  21H9G  JMD  APY  NODE  NU98ES 

2  »  10  X  »  1 21 6  ) 

<0G  5  1 = 1 » N 8 N 
I1=N0N(I> 

2FGCI )=7P(  II  ) 

X  ?  J (  I )=XPC  I  I  ) 

5  y=G( I )=YP(  II  ) 

Th=NPGWS-l 
CtTA=2./TP 
TP=NCQL-1 
OSI=  2. /TP 

CCi  2I  =  !#NP0WS 
T  F  *  I  ■  I 

E  IA  =1  •”TP*OE  TA 
CCI 2JS 1»NC0L 
T  f-J“l 

SI  =-l  .«T«*OSI 

1FCN  JK.E'K  M  00  TO  1012 

NCI  )  =  (  l.-<*  1  >*t  I. -FT  A)*<-1  >..♦?.•<  SI*  SI  4ETA  «ETA))/32. 
,Y(2)=S.*Cl.-ETA)*<l.-SI*Sl)*(  l.-3.*SI  >/32. 

N< J)  =  S.*(1 .-ETA >*C1.-SI*SI >*(  !.♦ 3.* SI  )/3  2. 

N(4  )  =  (  i.«S  I  )*C1.-ETAI*C-10.«9.*CSI  *SI  ♦  ETA*ETA)  )/32. 
N(5>=9.*U.#S1)*C1.-ETA*EIA)*C1.-3.*ETA)/32. 

N  (6)  =S.*C1  •*  SI  >*Cl#-tTA*ETA  )*  <1.«  3.*ETA  )  /  32. 

N<7  1 )*C I •«£TA)*C -10.«9 .*(  SI*SI«CT A •FT  A) >/32. 

NC8  .«ET A)*CI .-SI*SI )•< !•«  3.*S! )/32. 

N{9)=9.*Cl.«ETA)*(l,-S!*SI)*Cl.-3.*SI)/32. 

N(10)  =  U.-Sl)*(l.*erA)*(-10.«9,*CSl*Sl*ETA*£TA))/32. 
N(U)*9.*(1.-5I)*U.-ETA*ETA)*(1.4  3.*ETA)/32. 

*1  (12  >=y.*{  l.-SI  )  *(  l.-"TA*ETA  )*Ci.-3.*ETA)/32. 

GC  TO  1013 
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1 C 1 2  CONTINUE 

K<1  >  =  -'i.?5*(l.-Sl  >*Q.-ET4>*CSI*£l44l  .) 
M (2 ) -0.50*  n .-SI **2)*( l.-ET A) 

N(T)=C.25* (l.* SI  >*Cl.-tTA>*< SI-ETA-1.  ) 
N(4)- j  .50*  (1.4SI)*(t.-GTA**2) 

NC5) -Q .25* (1 .«S I)*  (J. 4  ETA)* (51 4  ETA-1. > 
N(5) -C. 50* (1.-51 **2)*{ I.* ETA) 
\(7)=C.25*(l.-Sl )*(1.«ETA)*(EIA-51-1. > 
\C  >=C.50*  (1  .-S!  )*  <1.  -  ETA**?) 

1  Cl  3  CONTINUE 

Z  C  ( I  .  J  )=  0. 

*C(!  * j)=0. 

red .  j>  =  o.o 

uC  l?  K=  l#  N-3N 

2  C  <  I  »J)  =  7C(I.J)4ZRG(K)*N(K) 

XC(!»J)  =  XCU.JMXRGCK>*N(K> 

12  YC(I . J >=YC (I »J >♦ YRG(K )*N(K> 

1M  =  1 
K  S  1  =  1 
KN2=NR  IwS 
KS2=NC3L 
CCS )  1-1.4 
NfiT=JT(N*3.I  ) 

If (Nfi  T  .EQ. Q.QR .NRT.GT.NHG)  GO  TO  50 
CC  36  J  =  t »  4 

56  IF( JTtNHT.  JJ.tQ.NRG)  N^TS-J 
K=NCOL 

tfd  .E0.2.0P.I.E0.4)  K=NFQWS 
Jl=l 

JK=rCOMPfI »  NRTS  ) 

If(JK.ra.-t)  JL  =  K 

CC  44  J  =  1.K 

GC  T (J  (45.46.47.46).! 

4  5  NN(N»OWS»J  )= NNRB ( NR T. NR T $. JL > 
KA2=NR3W5-l 
GO  TO  44 

46  NN(J*NC3L)=NNHQ(NRT .NRTS.JL) 

*S?=NCOL-t 

t  J  TO  4  4 

47  NM]  .  J  )=NMRB(NPT.NRTS.  JL  ) 

KM  =? 

GG  tq  uu 

43  M(  J.l  )=NNR=>(NH  T.NR  TS.  JL  > 

K  SI  =  2 

44  jL  =  JL 4  JK 
50  CCNUHUE 

IF  (KN1.GT.KN2)  GO  TJ  1'5 
IF  (KSl.GT .KS2)  GO  TO  105 
CC  10  I=KNl»KN? 

II=( I/?)*? 

OC  i )  J=KS 1 »*$? 


1 


r—  La 


J„=( J/2)*2 

IF  C  I. £3.  II.  AND.  J.fQ.  JJ)  GO  TO  t  C  01 
1 

MNU#  J)=NJ 
SC  TO  10 
100  1  NN(I»J)=) 
to  CCNTINUE 

00  4?  I  =1 »  NCOS. 

NNh’UNRG*!  *1  )  =  NN<NRQWS.I  > 

4 ?.  \NR3(NRG»3*I  )=NN(1»I) 

3G  4?  l=l/NRQ*S 
.NM-5<NRG»2*I  )=NN  ( I  *NCQL  ) 

43  NNR3CNRG»4  »!  )=NNCI*1  ) 

00  MO  I=l#NRO*S 
,0C  UC  J  =  1  » NCUL 
IF(NEfi-GT.NN(l»  J  )  )  GO  TO  210 
N£G=NN(I *  J  ) 
d  1C  CONTINUE 
IC5  CONTINUE 
K  =  1 

00  54  1=1- NFCWS 

CC  54  J= 1 / NCQL 

IFCNNC 1#J> .EQ.O)  GO  TO  54 

X  E ( K  )  =  XC( I »  J  ) 

YEH)*VC<I  *J) 

2E(K>  S7CC1  *  j) 

NECK  )  =  NN(  If  J  ) 

K  -  K+  1 

54  CONTINUE 
L=NP0*S-1 
00  151  I =1 *L»? 

00  151  *>  =  3-NC3L#2 

NFC1 )»(NC3L«<NCOL«l 1/2)4 J-? 

iF(  l ) =NR(l  )« l 
NR(J)=NR(2  )4  1 

NF(4)  =  NC0t  *Cl4l/2)4((NC'JL4l)/2>*(  I/2)4j-CJ/2> 

NFC5  )  =  <NC0L«  CNCOUl  >/?>*<I/2M  J 

Nfi(6)  =NH(5  )-l 

NS(M-NR<i)“l 

NFC5 )=NR(4 )-l 

NEL-NELU 

J1=N«<1) 

j2  =  NR( i  ) 

J3=NR  t  3  > 
o4  =  ‘jR(4  ) 

J5  =NR  (  5  > 

J6  =  NR<  6) 

J7  =  NR(7  ) 

J2 -»jr<t<3 ) 

L8U  )*1A35(N£1  J1  )-NE<  J?  )  )U 
LEC)s!43S(NE(  J2)-NE(  J3)H1 


L£ ( 3  )= I  43S  <NE<J3)“NEC  J4>  )-*l 
L  =  (4  )  =  1A3S  (NE(  J4  )-N£(J3i)4l 
L6(3  )  =  1  V3S<N£(  J5  >-NtA J6)  )4i 
L3(6)=IAJS(N£(  J6)-NEC  J7)  )4l 
L6t  /  )  =  1A0S(\'E<  J?  )-NEC  J«l  )4l 
L  E(? )  =  iAJSt NECJP  )-NECJl ) )4l 
1C  ?)?  IK=l,d 

If  (L  K  IX  ).  LE  .N3W  I  GO  T G  207 
NtW=Ld(IX> 

MEud  »«  =  NEL 
2 C?  CCNYINUE 

*>CNE(Jl  )>=XE(Jl  ) 

XKMECJ2))  =  X£<  J  2  ) 

<X(Mt(J3))=XE(J3) 

XXC*£C J 4  >)  =  X£( J4  ) 

*><M£(  J5))  =  X£( J  5  ) 

XXCNEC J6>»=X£(J6) 

XXt*E( J7  ))=XEl J7  ) 

>>CX“(  J6)>  =  X£(  J6  J 
YY<NE< J1 >> =YE( J1  ) 

YU'iZ(  J2>)=  Y£(J2) 

YY<NE<  Ji>)  =Y£{  J3  l 
YY(NE(J4)) =Y£< J4> 

YYU£(j5>>=y£CJ5> 

Y  Y(N£( J6 )) =Y£( J6  ) 

YY('JE(J7))=y£(J7) 

YY(.'|£(Jd))=YE(J8) 

22<M£{Jl )>=ZE(J!  ) 

2/<N£( J2>>  =  ZEC J2) 

ZZOI"(j3))  =  ZZ<J3) 

Z2(K~( J4  )  J  =  ZE< J4  > 

*/<MEl  J5)> =  <£( J5  ) 

ZZ(HE<JbM=ZE<  J6) 
lUHZi  J7  ))  =  Z£<  JT  ) 

Zif.NEl  Jd  ))  =  ZC<  Ja  ) 

N00(M£L#1)  =  N£(J1  ) 

NC9(N£L»2 )  =N  EC  J?  ) 

N CDC  NEL* 3 ) =  \E( J  3  ) 

NC0(N£l*4> =M£( J4  ) 

NCQCM£L»3>  =  N  EC  J5  ) 

NC0(NEi_»6)=N£CJ6  ) 

MC0CNEL.7) =NE< J7  ) 

NC0CNEL»aY=NECJ8  ) 

2  7c  CCrfMNUE 
151  CC.NTINLE 
16  CCNUNUE 

i 'SI  TEC»VK  )  MEL  *NEQ 
iGC  C0P*ATC2I5) 

CC  22  C  1  =  1 ,N£Q 

*MTECl*5'l>  XX{1)#YY(1)»ZZ(I)»U»V» 
22 C  CONTINUE 


END 

FILMED 
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